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ABSTRACT 


This  thesis  examines  the  vulnerability  of  the  Advanced  Encryption 
Standard  (AES)  to  algebraic  attacks.  It  will  explore  how  strong  the  Rijndael 
algorithm  must  be  in  order  to  secure  important  federal  information. 

There  are  several  algebraic  methods  of  attack  that  can  be  used  to  break  a 
specific  cipher,  such  as  Buchburger’s  and  Faugere’s  F4  and  F5  methods.  The 
method  to  be  used  and  evaluated  in  this  thesis  is  the  Multiple  Right  Hand  Sides 
(MRHS)  Linear  Equations.  MRHS  is  a  new  method  that  allows  computations  to 
be  more  efficient  and  the  equations  to  be  more  compact  in  comparison  with  the 
previously  referred  methods. 

Because  of  the  high  complexity  of  the  Rijndael  algorithm,  the  purpose  of 
this  thesis  is  to  investigate  the  results  of  an  MRHS  attack  in  a  small-scale  variant 
of  the  AES,  since  it  is  impossible  to  break  the  actual  algorithm  by  using  only  the 
existent  knowledge.  Instead  of  the  original  ten  rounds  of  AES  algorithm,  variants 
of  up  to  four  rounds  were  used. 

Simple  examples  of  deciphering  some  ciphertexts  are  presented  for 
different  variants  of  the  AES,  and  the  new  attack  method  of  MRHS  linear 
equations  is  compared  with  the  other  older  methods. 

This  method  is  more  effective  timewise  than  the  other  older  methods,  but, 
in  some  cases,  some  systems  cannot  be  uniquely  solved. 
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EXECUTIVE  SUMMARY 


This  thesis  compares  different  algebraic  methods  for  solving  and  breaking 
the  Advanced  Encryption  Standard  (AES)  algorithm.  However,  emphasis  is  given 
to  the  Multiple  Right  Hand  Side  (MRHS)  Linear  Equation  method,  which  is  a  new 
algebraic  method  used  for  attacking  ciphers  with  specific  algebraic  structures. 
Such  a  cipher  is  the  Rijndael  algorithm,  which  was  adopted  by  the  U.S. 
government  as  the  most  efficient  among  others  in  order  to  secure  federal 
information. 

In  cryptography,  it  is  considered  that  a  cipher  is  broken  when  a  method 
other  than  the  brute  force  attack  can  reduce  the  complexity  of  the  algorithm  or 
completely  decipher  it.  The  “brute  force  attack”  method  is  a  very  simple  method 
in  which  we  try  one  by  one  all  the  possible  keys  in  order  to  break  the  cipher.  In 
particular,  due  to  the  complexity  of  the  AES  algorithm,  it  becomes  very  difficult  to 
break  it  with  this  method.  Therefore,  AES  generated  interest  to  find  a  specific 
algorithm  to  break  this  cipher.  Some  of  the  older  algebraic  methods,  such  as 
Buchberger’s  algorithm  and  Faugere’s  F4  and  F5  algorithms,  which  are  briefly 
presented  in  Chapter  II,  cannot  break  this  cipher.  This  work  considers  a  new 
algebraic  method  (MRHS)  that  can  break  small  variants  of  the  AES.  Even  the 
breaking  of  a  small  variant  of  the  algorithm  constitutes  a  success,  since  it  opens 
a  new  horizon  in  the  cryptanalysis  field. 

A  comparison  with  the  older  algebraic  attack  methods  that  are  presented 
in  [9]  shows  that  the  MRHS  method  is  faster  than  the  others.  In  Chapter  II,  there 
is  a  quick  overview  for  some  of  these  methods  and  a  more  detailed  analysis  of 
MRHS  method  with  some  algebraic  examples. 

The  overall  concept  of  this  thesis  is  to  create  a  program  that  will  break  a 
small  variant  of  the  Rijndael  algorithm.  This  work  is  based  on  the  codes  that 
Professors  Havard  Raddum  and  Igor  Semaev  (of  University  of  Bergen,  Norway) 
provided  to  us.  These  codes  (in  the  “C”  language)  are  the  application  of  the 

xiii 


MRHS  method  to  the  AES  algorithm.  In  this  thesis,  we  created  new  programs, 
which  construct  the  equation  systems  of  the  AES  algorithm,  that  we  tried  to  solve 
by  using  the  codes  of  Professors  Raddum  and  Semaev.  Our  computational 
experiments  examined  many  cases  never  considered  in  previous  work. 

This  thesis  demonstrates  that,  for  a  small  variant  of  the  AES,  this  method 
is  very  efficient  (few  high-level  operations).  In  the  third  chapter,  some  examples 
are  contained  that  show  how  this  algorithm  works.  In  some  of  the  computational 
experiments,  a  totally  new  and  significant  result  appeared  for  the  cryptanalysis 
field;  sometimes  a  small  AES  variant  can  have  multiple  decryption  keys  for  a 
specific  encryption  key.  This  result  is  surprising  and  is  not  mentioned  anywhere 
by  Professors  Raddum  and  Semaev. 

In  particular,  the  first  example  shows  a  case  in  which  the  MRHS  method 
finds  multiple  solutions  to  the  system  (it  fails  to  solve  the  system — break  the 
cipher)  and  in  the  second  one,  all  the  MRHS  linear  equations  are  transformed 
into  ordinary  linear  equations  whose  total  number  is  the  number  of  variables  in 
our  system. 
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I.  INTRODUCTION 


A.  HISTORY 

In  1976,  the  United  States  government  as  an  official  Federal  Information 
Processing  Standard  (FIPS)  adopted  the  Data  Encryption  Standard  (DES). 
Today,  this  standard  is  considered  insecure,  and  the  necessity  for  a  new  more 
secure  and  efficient  standard  has  emerged.  Therefore,  in  1997,  the  National 
Institute  of  Standards  and  Technology  published  a  call  for  the  replacement  of 
DES.  Among  the  requirements  this  new  system  had  to  meet  were:  the  new 
algorithm  should  be  able  to  allow  key  sizes  of  128,  192  and  256  bits;  it  should 
operate  on  blocks  of  128  input  bits,  and  it  should  work  in  a  variety  of  different 
hardware.  For  example,  8-bit  processors  that  could  be  used  in  smart  cards  and 
the  32-bit  architecture  commonly  used  in  personal  computers  [1], 

For  the  adoption  of  the  new  algorithm,  five  finalists  were  chosen:  MARS 
(from  IBM),  RC6  (from  RSA  Laboratories),  Rijndael  (from  Joan  Daemen  and 
Vincent  Rijmen),  Serpent  (from  Ross  Anderson,  Eli  Biham  and  Lars  Knudsen) 
and  Twofish  (from  Bruce  Schneier,  John  Kelsey,  Doug  Whiting,  David  Wagner, 
Chris  Hall  and  Niels  Ferguson).  Finally,  in  2001,  the  Rijndael  algorithm  [2]  was 
adopted  as  the  Advanced  Encryption  Standard  (AES)  [3], 

B.  RELATED  WORK 

This  new  algorithm  was  of  great  interest  for  the  area  of  cryptanalysis.  The 
cryptanalists  started  to  seek  ways  to  attack  this  algorithm,  which  might  lead  to 
breaking  the  algorithm  (finding  the  key  more  efficiently  than  a  brute-force  attack). 
Since  Rijndael’s  algorithm  is  defined  using  algebraic  operations  in  finite  fields, 
cryptanalists  investigated  algebraic  methods  in  order  to  break  it.  Therefore,  they 
started  to  try  the  already  known  algebraic  techniques  for  solving  systems  of 
polynomial  equations.  Some  of  them  were  the  Buchburger’s  algorithm  and  the  F4 
and  F5  algorithms  [4]  that  are  summarized  in  Chapter  II  with  some  simple 
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examples.  Also,  they  tried  to  describe  the  new  algorithm  by  using  some  different 
algebraic  representations  over  the  Galois  finite  fields  [4]  that  are  discussed  in  the 
same  chapter. 

Recently,  a  new  algebraic  method  was  developed  that  used  Multiple  Right 
Hand  Sides  (MRHS)  Linear  Equations  [5],  The  main  difference  with  the  other 
algebraic  methods  is  that  the  equations,  which  describe  the  algorithm,  are  not 
expressed  anymore  as  multivariate  polynomial  equations,  but  instead  they  are 
presented  as  a  system  of  linear  equations,  each  having  a  set  of  multiple  right 
hand  sides. 

The  Norwegian  Professors  Havard  Raddum  and  Igor  Semaev  first 
developed  this  method.  After  a  comparison  between  the  Buchberger’s  F4 
algorithms  and  the  MRHS  method,  they  concluded  that  the  last  one  is  much 
faster. 

In  this  thesis,  after  we  borrowed  the  codes  [6]  from  the  Norwegian  team, 
we  managed  to  simulate  the  AES  algorithm  in  the  C  language  and  to  break 
small-scale  variants  of  the  Rijndael’s  algorithm.  After  the  application  of  the 
above-mentioned  codes,  many  computational  experiments  were  executed  to 
explore  how  the  method  would  work  in  many  different  cases.  Though  the 
previous  work  of  Raddum  and  Semaev  [5]  only  considered  AES  variants  using 
the  8-bit  field,  we  also  explored  variants  based  on  4-bit  and  2-bit  fields.  From 
these  experiments,  we  obtained  some  very  surprising  new  results  about  the 
decryption  (breaking)  of  variants  of  the  Rijndael  algorithm.  In  particular,  one  of 
the  them  is  that  there  are  some  cases  that  a  ciphertext,  which  was  encrypted  by 
a  specific  key,  could  be  broken  by  using  a  definite  number  of  different  keys.  This 
is  explained  in  detail  in  Chapter  III. 
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C.  INTRODUCTION  TO  AES 

Although  the  Rijndael  algorithm  was  adopted  as  the  official  algorithm  for 
AES,  it  is  important  to  note  that  the  algorithm,  which  is  used  in  the  AES,  is  a 
limited  version  of  the  Rijndael’s  algorithm.  In  particular,  Rijndael  is  a  block  cipher 
with  both  a  variable  block  length  and  a  variable  key  length  [2],  While  Rijndael’s 
algorithm  can  be  used  for  block  and  key  lengths  in  any  multiple  of  32  bits,  with  a 
minimum  of  128  bits  and  a  maximum  of  256  bits,  AES  uses  constant  block  length 
of  128  bits  and  key  lengths  of  128,  192  and  256  bits  only.  Also,  this  algorithm  is 
applicable  with  a  different  number  of  rounds  (Nr),  depending  on  the  number  of 
columns  of  the  cipher  key  (Nk)  and  the  number  of  columns  of  the  rectangular 
array  (Nb)  called  state,  which  is  actually  the  intermediate  cipher  result  of  the 
algorithm.  In  Table  1,  the  number  of  rounds  are  shown  as  a  function  of  the  block 
and  key  length. 


Nk 

Nb 

4 

5 

6 

7 

8 

4 

10 

11 

12 

13 

14 

5 

11 

11 

12 

13 

14 

6 

12 

12 

12 

13 

14 

7 

13 

13 

13 

13 

14 

8 

14 

14 

14 

14 

14 

Table  1 .  Number  of  rounds,  Nr,  as  a  function  of  Nb  (block  length/32)  and  Nk  (key 

length/32). 

Since,  in  this  thesis,  only  the  AES  algorithm  will  be  examined,  it  is  good  to 
analyze  the  subcases  of  the  Rijndael’s  algorithm  that  are  used  in  the  AES.  For 
the  encryption  process,  the  input  is  the  plaintext  block  and  the  initial  key  and  the 
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output  is  the  ciphertext  block.  In  the  decryption  process,  the  input  is  the 
ciphertext  block  and  the  output  the  plaintext  block. 

With  a  128-bit  key  length,  the  algorithm  applies  10  rounds.  (A  192-bit  key 
length  contains  12  rounds  and  a  256-bit  length  contains  14  rounds).  In  each 
round,  there  are  four  transformations  (linear  and  non-linear)  that  are  also  called 
layers.  Each  round  has  also  a  round  key,  derived  from  the  original  key  (input 
key).  The  round  transformation  and  its  steps  generate  intermediate  data  called 
states.  A  state  can  be  considered  as  a  rectangular  array  of  bytes  with  four  rows 
and  a  number  of  columns  (Nb)  that  depend  on  the  size  of  the  key  length.  Here, 
we  consider  that  the  key  length  is  128  bits,  where  the  key  is  arranged  in  a  4x4 
matrix  such  that  each  element  is  a  byte.  The  four  transformations  are  the 
ByteSub  transformation,  the  ShiftRow  transformation,  the  MixColumn 
transformation  and  the  AddRoundKey  transformation.  These  four  transformations 
compose  a  round.  The  four  transformations  are  discussed  next. 

D.  STRUCTURE  OF  THE  AES  ALGORITHM 

The  AES  algorithm  can  be  separated  in  two  stages — the  encryption  and 
the  decryption  process.  The  algorithm  for  the  encryption  process  includes  four 
transformations,  as  it  is  shown  below.  The  algorithm  for  the  decryption  process 
consists  of  the  inverses  of  the  above-mentioned  transformations  in  the  reverse 
order. 


1.  Encryption  Process 

The  sequence  of  the  four  transformations  mentioned  above  is  the 
following:  ByteSub,  ShiftRow,  MixColumn  and  AddRoundKey  transformations. 
We  next  present  these  transformations  in  order. 

a.  The  ByteSub  Transformation  (BS) 

This  is  the  only  non-linear  part  of  the  algorithm  and  assures 
resistance  to  differential  and  linear  cryptanalysis  attacks  [2],  This  transformation 

4 


consists  of  an  S-box,  which  is  applied  to  each  byte  element  of  the  state  (16-byte 
block)  independently  and  has  three  different  steps:  inversion,  a  Galois  Field  (GF) 
linear  mapping,  and  S-Box  constant  ,as  it  is  shown  in  Figure  1 . 
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Figure  1 .  The  ByteSub  step,  the  first  stage  in  a  round  of  AES  (From  [7]). 

(1 )  Inversion.  In  this  operation  of  the  S-box,  the  inverse  is 
computed  in  the  8-bit  Galois  Field,  GF(28).  The  byte  00000000  has  no  inverse 
and  00000000  is  used  in  place  of  its  inverse.  Assume  that  the  first  byte  is 
jc7x6x5x4jc3x2jcjX0  .  The  byte,  which  comes  up  from  the  inversion,  will  be 

v7  v6  v5 y4y3y2y1  y0  >  which  represents  an  eight  element  column  vector,  with  the 
rightmost  binary  bit  y0  in  the  top  position.  This  operation  provides  resistance 
against  the  linear  and  differential  cryptanalysis  attacks  [1], 

(2)  GF — Linear  Mapping.  At  this  point,  the  y  vector  is 
multiplied  by  a  constant  matrix,  and  the  column  vector  (0,1, 1,0, 0,0, 1,1)  is  added, 
yielding  a  vector  z7z6z5z4z3z2zjz0 :  (Note:  Galois  addition  is  equivalent  to  bitwise 
XOR  in  any  finite  field  of  even  size). 
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(3)  S-box  Table.  The  byte  z  is  the  input  to  the  S-box 
table.  Consider  an  input  byte  abcdefgh.  The  S-box  is  a  16x16  matrix.  We  look  for 
the  entry  at  abed  row  and  efgh  column  (rows  and  columns  are  numbered  from  0 
to  15).  The  intersection  of  these  two  entries,  transformed  into  a  binary  number,  is 
the  output  from  the  S-box.  For  example,  if  the  figure  10101010  is  the  input  byte, 
the  first  four  bits,  1010,  represent  the  decimal  10.  Therefore,  one  enters  at  the 
eleventh  row  and  eleventh  column.  The  intersection  is  172,  as  it  is  shown  in 
Figure  2.  This  is  converted  into  binary,  which  is  10101100.  This  is  shown  as  ‘ac,’ 
which  is  hexadecimal  for  10101100  in  Table  2.  That  number  is  the  output  of  the 
S-box. 
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Figure  2.  S-Box  in  the  encryption  process  (From  [3]). 
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In  order  to  make  clear  the  hexadecimal  notation,  Table  2 
shows  the  correspondence  between  decimal,  binary  and  hexadecimal  notations. 


Decimal 

Binary 

Hexadecimal 

0 

0000 

0 

1 

0001 

1 

2 

0010 

2 

3 

0011 

3 

4 

0100 

4 

5 

0101 

5 

6 

0110 

6 

7 

0111 

7 

8 

1000 

8 

9 

1001 

9 

10 

1010 

A 

11 

1011 

B 

12 

1100 

C 

13 

1101 

D 

14 

1110 

E 

15 

1111 

F 

Table  2.  Correspondence  between  decimal,  binary  and  hexadecimal  notations. 

b.  The  ShiftRows  Transformation  (SR) 

In  this  transformation,  which  is  linear,  the  rows  of  each  state  are 
cyclically  shifted  to  the  left,  with  each  row  shifted  a  different  amount,  as  it 
appears  in  Figure  3.  This  provides  resistance  against  truncated  differential  and 
saturation  attacks  [2],  For  example,  row  zero  is  shifted  by  Co  bytes,  row  1  is 
shifted  by  C-i  bytes,  row  2  is  shifted  by  C2  bytes  in  such  a  way  that  the  byte  at 
position  j  in  row  i  moves  to  position  (j-Cj)  mod  Nt>.  Particularly  in  AES  C0  =  0 , 

Cj  =1,  C2  =2  and  C3  =  3.  The  output  of  this  transformation  is  a  matrix  with  the 
same  dimensions  as  the  input  matrix. 
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Figure  3.  The  ShiftRows  step,  the  second  stage  in  a  round  of  AES  (From 

[7]). 


c.  The  MixColumns  Transformation  (MC) 

This  transformation  operates  on  each  4-byte  column  separately  and 
is  omitted  in  the  last  round.  It  is  also  a  linear  transformation,  which  has  diffusion 
power.  The  columns  of  the  state  are  considered  as  polynomials  over  GF(28), 
which  are  multiplied  by  a  fixed  polynomial  c(x)  modulo  (x4+1).  This  polynomial  is: 

c(x)  =  03  ■  x3  +  0 1  •  x2  +  0 1  ■  x  +  02  [2] 

and  this  multiplication  can  be  presented  as  a  matrix  multiplication  [2]  (in 
hexadecimal  notation),  as  it  appears  below: 


K 

"02 

03 

01 

or 

«o 

k 

01 

02 

03 

01 

X 

d\ 

b2 

01 

01 

02 

03 

a2 

b3  _ 

03 

01 

01 

02_ 

a3 

Figure  4  depicts  the  MixColumns  operation  on  the  columns  of  the 

state. 
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Figure  4.  The  MixColumns  step,  the  third  stage  in  a  round  of  AES  (From  [7]). 


d.  The  AddRoundKey  Transformation 


In  this  transformation,  a  key,  consisting  of  128  bits,  which  are 
arranged  in  a  4x4  byte  matrix,  is  added  to  the  output  of  the  MixColumn 
transformation.  A  different  round  key  is  added  to  the  state  at  the  end  of  each 
round. 
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Figure  5.  The  AddRoundKey  step,  the  fourth  stage  in  a  round  of  AES  (From 

[7]). 


This  key  is  derived  recursively  from  the  original  key  as  follows.  We 
will  see  this  procedure  in  five  steps. 

Stepl:  Label  the  first  four  columns  of  the  original  key  W(0),  W(1), 

W(2),  W(3). 
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Since  the  whole  algorithm  consists  of  10  rounds,  40  more  columns 
are  required,  four  for  each  round.  Let  i  symbolizes  the  number  of  column  in  the 
different  round  keys  as  the  columns  are  derived  from  the  key  schedule 
(4 < z  <  43 ). 

Step  2:  If  the  number  of  the  new  column  is  a  multiple  of  four, 

then... 


Step  3:  It  is  W(i)  =  W(i-4)®T{w(i-\)),  where  T(w(i- 1))  is  the 

transformation  of  W(i- 1)  obtained  as  follows.  Assume  that  the  elements  of  the 

column  W(i- 1)  are  a,b,c,d.  These  are  shifted  cyclically  to  obtain  b,c,d,a,  and 

then,  at  this  point,  each  of  these  bytes  is  substituted  with  its  corresponding  byte 
from  the  S-box  of  the  ByteSub  transformation,  as  it  is  explained  above.  So,  four 
other  bytes  result,  i.e. ,  e,f,g,h. 

Step  4:  Finally,  the  round  constant 

r  (/)  =  000000 10('  4)74 


is  computed  in  GF(28).  So,  the  T(w(i- 1))  is  the  column  vector  (e@r(i),f,g,h) . 


Step  5:  If  i  is  not  a  multiple  of  four,  then 

W{i)  =  W(i-A)®W(i-l). 

2.  Decryption  Process 

The  decryption  process  consists  of  the  inverses  of  the  four  encryption 
steps:  InvByteSub,  InvShiftRows,  InvMixColumns  and  invAddRoundkey. 
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a.  The  InvByteSub  Transformation 

This  transformation  consists  of  the  inverse  S-box.  In  essence,  in 
this  step,  the  inverse  transformation  of  the  equation  that  was  made  in  the 
ByteSub  transformation  is  performed.  For  the  linear  mapping,  one  takes: 
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The  inverse  of  the  S-box  appears  in  Figure  6: 
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Figure  6.  S-Box  in  the  decryption  process  (From  [3]). 


11 


b.  The  InvShiftRows  Transformation 

In  this  transformation,  the  opposite  shifting  operation  is  applied. 
Therefore,  the  rows  are  shifted  to  the  right  instead  of  to  the  left,  which  takes 
place  at  the  ShiftRows  transformation. 

c.  The  InvMixColumns  Transformation 

In  this  transformation,  every  column  is  multiplied  by  the  inverse 
polynomial  of  c(x)  (mod  x4+1 )  which  is: 

d  (x)  =  OB  ■  x3  +  0 D  •  x2  +  09  •  x  +  0 E 

The  inverse  matrix  multiplication  of  the  equation,  which  was  used  in  the 
MixColumn  transformation,  is: 
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This  transformation  is  omitted  in  the  last  round. 

d.  The  InvAddRoundKey  Transformation 

This  transformation  applies  the  keys  that  were  used  in  the 
encryption  process  in  the  reverse  order. 

E.  THESIS  OBJECTIVE 

A  common  method  to  break  a  cipher  is  called  “the  brute  force  attack.”  In 
this  method,  one  tries  to  break  a  cipher  using  one  or  more  ciphertexts,  where  the 
corresponding  plaintexts  are  known,  by  trying  to  decrypt  them  using  all  the 
possible  keys  for  the  cipher.  When  the  decryption  yields  the  known  plaintext,  we 
have  found  the  key  and  the  cipher  is  broken.  (A  variation  on  this  method,  used 
when  the  plaintexts  are  not  known,  stops  when  a  meaninful  plaintext  is  found). 
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However,  as  ciphers  become  more  and  more  complicated,  this  method 
turned  out  to  be  useless.  For  example,  in  the  case  of  the  Rijndael’s  algorithm, 
one  has  to  try  at  most  2128  possible  keys  to  break  the  cipher.  This  is  impractical, 
since  the  time  and  the  memory  that  are  necessary  for  that  exceed  the  limits  of  all 
existing  workstations. 

In  the  last  few  decades,  some  cryptanalysis  has  been  based  on  algebraic 
attacks.  Algebraic  attacks  have  the  advantage  that  a  cryptanalyst  does  not  need 
to  have  many  known  plaintexts  and  ciphertexts  in  order  to  create  an  equation 
that  can  describe  the  cipher,  an  event  that  happens  with  the  linear  or  differential 
cryptanalysis  where  one  has  to  have  many  pairs  of  plaintexts  and  ciphertexts  in 
order  to  be  able  to  describe  the  cipher. 

In  this  thesis,  a  specific  algebraic  attack  method — the  Multiple  Right  Hand 
Side  (MRHS)  Linear  Equations  method — is  examined  and  compared  with  other 
older  methods  (Buchburgers,  F4,  F5)  [4],  In  addition,  a  small-scale  variant  of  the 
AES  will  be  examined,  since,  as  it  will  be  explained  later,  the  number  of  the 
variables  for  the  whole  AES  algorithm  is  too  big  to  be  solved. 

Therefore,  in  Chapter  II  we  will  briefly  summarize  these  older  algebraic 
methods.  While  these  methods  can  represent  the  AES,  they  cannot  break  the 
cipher.  A  more  detailed  analysis  of  the  MRHS  Linear  Equations  method  will  also 
be  presented  with  a  very  simple  example  in  order  to  show  how  this  method 
works. 

In  Chapter  III,  we  will  present  the  computational  experiments,  where  we 
apply  the  MRHS  approach  to  small  variants  of  AES,  including  several  that  have 
never  been  attacked  before  by  this  method.  We  present  the  results  of  the 
application  of  this  method  with  two  representative  examples.  In  the  first  example, 
the  method  fails  to  solve  the  algebraic  system  that  is  created  with  the  MRHS 
Linear  Equations,  while  in  the  second  one  we  conclude  to  a  solution. 
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II.  ALGEBRAIC  EQUATIONS  FOR  AES 


A.  INTRODUCTION 

There  are  many  different  ways  to  describe  a  cipher.  However,  the 
complexity  of  modern  ciphers  requires  knowledge  of  their  algebraic 
representation  in  order  to  attack  and  break  them.  In  particular,  one  tries  to 
describe  a  cipher  by  finding  the  algebraic  properties  that  it  has  and,  after  that,  by 
creating  some  homomorphisms  or  isomorphisms  of  this  cipher.  These  new 
structures  are  very  helpful,  both  to  the  implementers  of  a  particular  cipher  who 
want  to  provide  further  protection  to  the  cipher  against  side-channel  attacks,  and 
to  the  cryptanalysts  who  try  to  analyze  it  or  even  to  break  it. 

B.  DIFFERENT  REPRESENTATIONS  OVER  F2  AND  F256 

The  meaning  of  F2  is  the  finite  field  with  order  (number  of  elements) 
pn  =  21 ,  where  p  is  always  a  prime  number  and  n  is  a  positive  integer. 
Respectively,  F256  is  the  finite  field  with  order  pn  =2X  =256  ,  which  is  also  notated 
GF(28). 

After  the  adoption  of  the  modern  algorithms,  one  of  which  is  the  Rijndael’s 
algorithm,  a  great  interest  was  created  for  the  wider  application  of  computational 
algebra  in  cryptography.  Therefore,  a  number  of  different  algebraic 
representations  were  developed. 

In  the  case  of  AES,  one  can  consider  a  number  of  different 
representations,  which  are  called  dual  ciphers  depending  on  the  properties  of 
their  representation  mappings.  Some  are  described  in  the  next  pages. 

1.  Isomorphic  Representations 

Definition:  Suppose  A  is  a  vector  space  over  a  field  F  with  a 
multiplication  operation  Ax  A  —>  A .  If  this  multiplication  operation  is  associative 

15 


and  is  a  bilinear  mapping  on  the  vector  space  A ,  then  A  is  an  (associative)  F  - 
algebra,  or,  more  simply,  an  algebra[ 4], 

In  these  ciphers,  the  mappings  of  the  state  and  key  spaces  are  algebra 
isomorphisms  of  the  AES  state  space  algebra,  where  algebra  is  defined  below. 
Therefore,  these  ciphers  are  isomorphic  to  the  AES. 

In  AES,  each  byte  can  be  considered  as  an  element  of  the  finite  field 
GF[2 in  terms  of  the  following  polynomial: 

m(x)  =  x8  +x4  +x3  +x  +  l  (Rijndael  polynomial) 

which  is  irreducible  in  GF( 2)[x] .  This  finite  field  can  be  constructed  in  many 
different  ways  from  the  chain  of  its  subfields,  as  it  appears  below: 

GF (2)  <=  GF (22 )  c=  GF (24 )  <z  GF (2s ) 

In  total,  there  are  30  +  (l  •  60)  +  (3  •  120)  +  (1  -  6  - 120)  =  1 170  different 
isomorphic  representations  of  the  AES  based  on  the  different  irreducible 
polynomials  of  the  GF[2 subfields  as  they  appear  in  Table  3  [4], 


Subfield 

Degree 

GF  (2) 

GF  (22) 

GF  (24) 

2 

1 

6 

120 

4 

3 

60 

- 

8 

30 

- 

- 

Table  3.  The  number  of  irreducible  polynomials  over  subfields  of  GF  (28). 


These  particular  representations  are  not  of  cryptanalytic  interest.  Rather, 
they  are  intended  to  improve  the  efficiency  of  hardware  implementation. 
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2.  Regular  Representations 


The  regular  representation  is  the  algebra  homomorphism  v:  A—>  Mn  (iT) 

that  maps  a  e  A  to  the  matrix  corresponding  to  the  linear  transformation  z  i-^  az , 
where  z  is  a  vector  over  F  of  length  n  [4], 

3.  Logarithmic  Representations 

Since  an  element  of  a  finite  field  can  be  represented  with  logarithmic 
functions  instead  of  vector  spaces,  an  element  of  the  AES  state  space  (excluding 

zero  bytes)  can  be  described  as  an  element  of  the  set  (Z255  )16 .  These  ciphers  are 

called  log  dual  ciphers  and  there  are  128  different  primitive  elements  in  F  ,  giving 
128  such  representations  [4],  More  details  of  how  to  specify  a  logarithmic 
representation  of  the  AES  are  given  in  [8], 

C.  ALGEBRAIC  SOLUTION  METHODS 

In  this  subsection,  some  algebraic  solution  methods  will  be  discussed.  The 
first  two  will  be  based  on  the  Grobner  basis  algorithms,  which  are  well-known 
methods  for  the  solution  of  multivariate  polynomial  equations.  These  two 
algorithms  have  been  unable  to  break  the  Rijndael  algorithm,  and  at  the  end  of 
this  chapter  there  will  be  a  comparison  between  these  algorithms  and  the  MRHS 
Linear  Equation  algorithm  in  terms  of  their  complexities  [4],  The  last  one,  which  is 
examined  and  analyzed  in  this  thesis  with  more  details,  is  a  new  approach  that 
seems  to  have  better  results  than  the  previous  two  based  on  the  time  that  is 
needed  in  order  to  break  the  AES  cipher. 

1.  Buchberger’s  Algorithm 

As  mentioned  above,  Buchberger’s  algorithm  solves  systems  of 
multivariate  polynomial  equations.  Consider  a  polynomial  ring  F[jcj,x2,...,jc„]  with 
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a  monomial  ordering.  Suppose  I  <  F[x15x2,...,xn]  is  an  ideal  of  this  ring  with  a 
basis  F  =  .  The  S-polynomial  of  any  pair  of  the  function  of  the  basis  is 

defined  as: 


where  Icm  is  the  least  common  multiple,  LM  is  the  Leading  Monomial,  and  LT  is 
the  Leading  Term.  What  one  tries  to  achieve  with  this  polynomial  S,  which 
belongs  to  the  ideal  I ,  is  to  cancel  the  leading  terms  from  any  pair  of  the 
polynomials  and,  with  the  help  of  the  next  theorem,  to  compute  a  Grobner  basis 
of  the  ideal  I  [4], 

Theorem:  Let  F[xpx2,...,x„]  be  a  polynomial  ring  with  a  monomial  ordering,  and 
let  I  be  an  ideal  of  ¥\xl,x2,...,xn\.  A  basis  G  =  for  the  ideal  I  is  a 

Grobner  basis  for  I  if  and  only  if  every  S-polynomial  s(fnfj)  of  pairs  of  distinct 
polynomials  fi,fJeG  has  remainder  0  upon  division  by  G  . 


Consider  an  example  of  how  to  solve  a  multivariate  polynomial  system  by 
using  this  algorithm. 

Example:  Consider  a  polynomial  ring  F[x,y]  with  multivariate  polynomials 
(with  two  variables)  over  the  complex  numbers  with  the  lexicographic  order 
(v-<x).  The  ideal  in  this  particular  example  is  generated  from  the  polynomials 
that  appear  below: 

f  =  x2y- 1  and  f2  =  xv2-x 


Here,  the  Grobner  basis  is  computed  by  using  the  Buchburger’s  algorithm.  Set 
G  —  {f ,  f i }  and  compute  the  S(f,f2)  -polynomial. 


2  2 

x  y  (  2 
-{x  y- 


2 

x  y 


i) 


2  2 

x  y 


xy 


18 


=  y  (x2y  ~ l)  - x(xy2  -  x)  =  x2  -  v 

The  leading  term  of  the  S(fl,f2)  is  x2.  The  leading  term  of  /j  and  f2  is  x2v  and 
xy2,  respectively.  So,  S(fl,f2)  cannot  be  reduced  by  these  two  polynomials.  The 
G  is  expanded,  so  that  it  becomes  G  =  {/1,/2,/3}  where  f3=x2-y  and  ’) 

and  S(f2,f3)  are  computed. 

S{fl,fi)  =  y2-\,  which  cannot  be  reduced  by  the  set  G  =  {fx,f2,f3} .  So  the  set 

G  =  {f,f2,f3,f4}  is  expanded,  where  f4=y2- 1. 

S(f2,f3)  =  -x2+y\  Its  leading  term  is  -x2  which  is  divisible  by  the  LT{f3).  The 
division  of  S(f2,fi)  by  /3  and  then  by  f4  results  in: 

S{f2Ji)  =  -x2  +  v3  =  -f3  +  ( v3  -  v)  =  -/3  +  yf4 

Now,  S(f\,f4)=  f3,  S(f2,f4)  =  0  and  S(f3,f4)=f3-yf4  are  also  computed.  Thus, 
one  concludes  that  all  these  S-polynomials  can  be  reduced  by  G.  So, 
G  =  {/1;/2,/3,/4}  is  a  Grobner  basis  of  the  ideal  /  .  The  reduced  Grobner  basis  is 

G  =  |x2  —  y, y2  —  1 J  ■  Finally,  the  following  system  of  equations  is  solved  in  order  to 
find  the  solution. 

x2  -  v  =  0 
v2  - 1  =  0 

The  complete  solution  of  this  system  is  {(1, 1) ,(—l,l),(z',—l) 1)} . 

2.  F4  and  F5  Algorithms 

These  two  algorithms  have  their  names  from  their  creator,  Jean  Charles 
Faugere,  and  compute  again  a  Grobner  basis.  They  can  be  viewed  as  an 
improved  method  of  Buchburger’s  algorithm.  They  are  based  on  the  same 

principles.  In  particular,  the  F4  algorithm  instead  of  polynomial  reduction  uses 
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matrix  reduction,  while  F5  uses  matrix  reduction,  but  each  of  these  generated 
matrices  is  of  full  rank.  This  method  is  illustrated  in  the  next  example  [4], 


Example:  Consider  the  polynomial  ring  R[x,y,z],  which  includes 

polynomials  with  three  variables  in  lexicographic  order.  We  want  to  reduce  the 
following  polynomials: 

f  =3x  yz-5xy  and  f2  =  5x2z2+3xv  +  l 

by  the  set  of  polynomials  {g^g.,}  ,  where 

gj  =xy-2z  and  g2  =  x2z-3yz. 

To  reduce  f\ ,  with  respect  to  {gpg2} ,  the  following  reduction  is  performed: 

f  =  3x  yz  -  5xy 

=  6  x1  z1  -5xy +  (3x2z)g, 

=  -5xy  + 1 8yz2  +  ( 3x2z )  ^  +  ( 6z )  g2 
=  18yz2  -10z  +  (3x2z)gj  +(6z)g2  -(5)gj 

So,  /  is  reduced  with  respect  to  {gpg2}  to  18yz2-10z. 

The  same  procedure  for  f2  is 

f2  =  5 x  z2  +  3 xv  + 1 

=  3xv  +  15yz2  +l  +  (5z)g2 
=  15yz2  +  6z  +  l  +  (5z)g2+(3)g1 

Thus  f2  is  reduced  to  15yz2+6z  +  l  with  respect  to  {gpgj. 
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The  idea  behind  the  F4  and  F5  algorithms  is  to  make  these  reductions  as  a 
matrix  reduction.  However,  the  reduction  for  both  f\  and  f2  requires  reduction 

only  with  respect  to  (Vz)^,  g, ,  and  ( z)g2 .  Therefore,  a  matrix  of  coefficients  is 

created. 

3  2  2  2  1 

x  yz  x  z  yz  xy  z  1 

fx  (3  0  0  -5  0  0 

f2  0  5  0  -3  0  1 

x2zgl  1  -2  0  0  0  0 

lg,  0  0  0  1  -2  0 

zg2  v0  1  -3  0  0  0 

The  reduction  steps  correspond  to  the  row  reduction  of  the  first  two  rows, 
which  represent  the  polynomials  f\ ,  f2  using  the  last  three  rows,  which  represent 

the  (x2z]gj,gj,(z)g2  polynomials.  So,  as  a  result  of  the  reduction,  the  following 

matrix  is  derived. 

x  yz  x  z  yz  xy  z  1 

0  0  18  0  -10  0 

0  0  15  0  6  1 

1  -2  0  0  0  0 

0  0  0  1  -2  0 

0  1  -3  0  0  0 

The  first  two  rows  of  the  reduced  matrix  give  the  reduction  of  / ,  f2  with  respect 
to  g1?g2 .  Note,  that  this  is  the  same  result  that  was  achieved  with  the 
Buchburger’s  reduction. 


A 

A 

x2zgl 

l£i 

Z&2 
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3.  Multiple  Right  Hand  Sides  (MRHS)  Linear  Equations  Algorithm 


A  different  approach  of  algebraic  attack  not  based  on  equations  of 
polynomials  will  be  examined.  This  presentation  will  utilize  a  special  type  of 
equation  known  as  Multiple  Right  Hand  Side  Linear  Equations. 

Definitions: 

1 .  A  set  of  small-scale  Variants  of  the  AES  denoted  by  SR(n,r,c,e)  are 
defined,  where 

•  n  is  the  number  of  rounds, 

•  r  is  the  number  of  rows  in  the  rectangular  arrangement  of  the  input, 

•  c  is  the  number  of  columns  in  the  rectangular  arrangement  of  the 
output, 

•  e  is  the  number  of  bits  in  a  “word.” 

Since  we  consider  that  e=8,  the  underlying  finite  field  is  GF(28)  and  all  the 
matrices  and  vectors  are  over  GF(2). 

2.  Let  X  be  a  set  of  Boolean  variables  represented  as  a  column 
vector.  An  equation  of  the  form 

AX  =  al,a2,...,as 

or 

S,:4X  =  [Z,] 

is  called  an  MRHS  system  of  linear  equations  if  A  is  a  matrix  of 
size  kxn  and  rank  k  and  ava2,...,as  are  column  vectors  of  length  k. 

3.  A  symbol  S=(X,L)  consists  of  an  ordered  set  of  variables  X=X(S) 
and  a  list  L=L(S)  of  Right  Hand  Sides  [9] 

a.  Agreeing  Procedure 

One  of  the  problems  with  applying  algebraic  attacks  in  a  cipher  is 
the  difficulty  of  presenting  practical  examples,  because  the  required  time  and 
memory  requirements  grow  beyond  the  limitations  of  a  typical  workstation. 
Therefore,  only  small-scale  variants  of  the  different  ciphers  were  performed. 
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After  the  completion  of  the  rounds  of  AES,  one  can  count  the 
number  of  variables  there  are  in  the  system  in  order  to  determine  the  algebraic 
equations  that  can  describe  the  system. 

For  example,  the  total  number  of  variables  for  the  AES,  assuming 
8-bit  words,  is  given  by  the  following  formula: 

8  rc  +  8  nr  +  8(»  - 1  )rc 

Thus,  for  a  complete  AES  the  number  of  variables  is  1600.  This  is  large,  if  one 
wants  to  solve  a  system  with  so  many  non-linear  polynomial  equations. 

Figure  7  shows  the  bytes  that  are  variables  in  a  small-scale  variant 

of  the  AES. 


Kp  Kt  Kj  K3 


Figure  7.  SR(3,4,4,8)  equations  variables  (From  [9]). 

The  variables  can  be  separated  in  two  different  categories.  The  first 
one  includes  all  the  bits  that  are  derived  after  each  S-box,  except  the  last  one 
that  is  considered  as  known,  since  the  ciphertext  is  considered  as  a  known  for 
this  algorithm.  The  second  one  is  the  16  bytes  of  the  initial  key  and  after  that, 
only  the  first  column  of  each  expansion  key,  since  the  first  column  is  derived  with 
the  contribution  of  the  S-box,  while  the  three  others  are  dependent  on  the  first 
column,  as  it  is  shown  in  Figure  7  and  described  in  the  AddRoundKey 
transformation  in  Chapter  I. 
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Consider  two  symbols: 

Sr.A,X  =  [L]  and  S^.AJ  =  [Lt] 

The  matrices  L,  are  of  size  /c,xs;.  These  symbols  are  derived  from 
the  bits  that  go  in  and  out  of  each  S-box,  which  are  expressed  as  linear 
combinations  of  the  variables. 

These  two  symbols  agree  if,  for  any  ax  e  L. ,  there  exists  an  a2  e  Z, 
such  that  the  linear  system 


4 

x  = 

a\ 

aj\ 

_a2  _ 

is  consistent,  and,  conversely,  for  any  a2  el.,  there  exists  an  ax  e  T  such  that 


the  linear  system  is  consistent.  The  steps  for  the  application  of  this  attack  are  as 
follows: 


Define  the  matrices: 


"4" 

\L~ 

"0  " 

A  = 

,T„  = 

1 

,T„  = 

aj\ 

9  u 

0  _ 

5  Ji 

LJ_ 

with  t  =  kt  +kj  rows. 


•  Choose  a  nonsingular  transform  matrix  U  of  size  txt  such  that  the 
product  UA  is  a  matrix  with  zeros  in  its  last  r  rows.  If  r=0,  then  the 
symbols  agree. 

•  If  r>0,  then  compute  the  matrices  UT..  and  UT.. .  Let  Pr  denotes 
the  set  of  UT.. -column  projections  to  the  last  r  coordinates.  If 
Pr&.  =  Pr.  the  symbols  agree. 

•  If  Pr  *Pr..  then  remove  the  columns  from  L.  whose  image  is  not 
found  in  Pr^.flPr,.,.  and,  similarly,  the  columns  from  L.  whose  image 
is  not  found  in  Pr.fl  Pr/(. .  One  concludes  with  two  new  symbols 
whose  equations  are  S\ :  AiX  =  [4]  and  S) :  AjX  =  [z\] . 


This  procedure  is  repeated  up  to  the  point  that  all  symbols  agree.  When  the 
agreeing  procedure  is  applied  in  a  pair  of  symbols  S,  and  Sy  and  the  procedure  is 
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continued  with  other  symbols  in  the  system,  for  example,  Sk ,  it  is  possible  that 

the  Si  and  Sk  will  disagree.  This  means  that  S,  and  Sj  may  disagree  again  after  S, 
and  Sk  agree.  So,  the  previous  agreements  have  to  be  run  through  again  and 
again. 


b.  Gluing  Procedure 

After  the  completion  of  the  agreeing  procedure,  one  may  conclude 
that  the  system  of  equations  does  not  have  a  solution.  Therefore,  a  way  to  start 
the  agreeing  procedure  again  must  be  found.  Here  is  where  the  gluing  procedure 
comes  into  play.  This  method  merges  two  symbols  into  one,  bringing  their  joint 
information  in  the  new  symbol.  The  agreeing  procedure  is  applied  again  with  the 
new  symbols,  and  is  repeated  until  a  unique  solution  is  reached.  The  following 
example  demonstrates  how  this  works.  Also,  this  step  increases  the  complexity 
of  the  attack. 


c.  Example 


1.  Agreeing  Procedure.  We  have  two  equations 

A,X  =  [L,]  and  A2X  =  [L2]  in  variables  X  =  {xvx2,x2,x4,x5} : 


*1 

*1 

“l 

1 

0 

0 

o' 

*2 

"l 

0 

0 

ll 

"0 

1 

0 

0  1" 

x2 

"0 

1 

0 

o' 

1 

0 

1 

0 

0 

x3 

= 

0 

1 

0 

0 

J 

0 

0 

1 

0  1 

x3 

= 

1 

1 

0 

0 

1 

0 

0 

1 

0 

*4 

0 

0 

1 

1 

0 

0 

0 

1  1 

*4 

1 

1 

0 

1 

_x5J 

_X5J 

By  following  the  previously  described  procedure,  the  matrix  A  is: 

"1  1  0  0  0" 

10  10  0 

10  0  10 

A  = 

0  10  0  1 
0  0  10  1 
0  0  0  1  1 
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A  non-singular  matrix  U  for  the  transformation  of  A  is: 

"1  0  0  0  0  0 

0  1  0  0  0  0 

0  0  1  0  0  0 

u  = 

0  0  0  1  0  0 

110  110 
10  110  1 

So, 

"11000 
10  10  0 
10  0  10 

UA  = 

0  10  0  1 
0  0  0  0  0 
0  0  0  0  0 

Then,  r=2. 
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Looking  at  the  last  two  rows  of  the  above  matrices,  one  defines: 

Pri2={(1 .1  ).(1 .0),(0,1 )},  Pr21={(1 ,1 ), (0,0), (0,1 )} 

Pri2nPr21={(1,1),(0,1)} 

One  can  observe  that  the  second  and  the  fourth  column  of  UTn  do  not  match 
with  any  column  of  UT1X.  Therefore,  these  columns  are  removed.  Similarly,  the 
second  and  the  third  column  of  UTn  should  be  removed  for  the  same  reason. 


In  addition,  two  new  symbols  result, 


*1 

*1 

"1 

1 

0 

0 

0" 

x2 

'1 

o' 

"0 

1 

0 

0 

f 

x2 

"0 

0" 

1 

0 

1 

0 

0 

x3 

= 

0 

0 

j 

0 

0 

1 

0 

1 

x3 

= 

1 

0 

1 

0 

0 

1 

0 

X4 

0 

1 

0 

0 

0 

1 

1 

X4 

1 

1 

_v 

_x5_ 

that  now  agree. 

2.  Gluing  Procedure.  Let  B  be  the  sub-matrix  of  UA  in  its 
last  t-r  nonzero  rows.  The  gluing  of  the  two  previously  agreed  symbols  is  BX=[L], 
where  each  column  of  L  is  the  sum  of  one  column  from  UT12  and  one  from  the 
UT2i  with  the  same  projection  in  its  last  r  coordinates,  reduced  to  the  first  t-r 
rows. 

For  the  above  example,  the  gluing  procedure  yields: 


"1 

1 

0 

0 

0" 

*1 

'1 

0" 

1 

0 

1 

0 

0 

x2 

0 

0 

1 

0 

0 

1 

0 

x3 

— 

0 

1 

0 

1 

0 

0 

1 

*4 

0 

0 

_x5J 

This  MRHS  system  of  linear  equations  contains  the  information  of  the  first  two 
different  symbols. 
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d.  From  MRHS  to  Linear  Equations 

In  order  to  derive  the  ordinary  linear  equations  (unique  right  hand 
side)  from  the  MRHS  linear  equations,  the  L  matrix  is  triangulated  with  a  row 
transformation.  An  upper-triangular  matrix  with  zeros  in  its  last  r-i>0  rows  is 
sought,  from  where  to  take  r-i  homogeneous  equations.  There  may  also  be  non- 
homogeneous  equations  with  ones  in  the  whole  row  of  L.  In  the  above  example, 
the  triangulation  of  L  results  in  the  following  symbol: 


"1 

1 

0 

0 

0" 

*1 

"1 

0" 

0 

1 

0 

1 

0 

x2 

1 

1 

1 

0 

1 

0 

0 

x3 

— 

0 

0 

0 

1 

0 

0 

1 

*4 

0 

0 

_v 

which  is  equivalent  with  the  three  following  linear  equations  and  the  two  initial 
MRHS  linear  equations  [5]: 

x2+x4  =  l 
Xj  +  x3  =  0 
x2  +  x5  =  0 

4.  Algorithms’  Complexities  and  Comparison 

This  section  compares  the  complexities  of  the  previously  described 
algorithms. 

The  Buchberger’s  algorithm  time  complexity  is  related  to  the  total  degree 
of  all  the  intermediate  polynomials  that  are  generated  by  the  algorithm.  In 
essence,  this  algorithm  can  have  double  exponential  complexity.  In  particular, 
considering  the  AES  equation  system  over  GA( 28),  the  complexity  of 
Buchburger’s  algorithm  is,  at  worst,  single  exponential  [4], 

The  F4  and  F5  algorithms  are  different  approaches  for  the  computation  of  a 
Grobner  basis.  Since  the  Buchberger’s  algorithm  involves  polynomial  reductions, 
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which  take  place  sequentially,  while  F4  and  F5  use  matrix  reductions,  it  seems 
logical  that  F4  and  F5  will  be  faster.  The  F5  algorithm  is  even  better  than  the  F4 
because  it  uses  only  full  rank  matrices. 

The  complexity  of  the  MRHS  algorithm  is  very  difficult  to  specify  since  the 
different  parameters  that  can  be  used,  such  as  the  number  of  rounds,  the 
number  of  rows  and  columns,  etc,  make  it  very  complicated.  However,  this 
algorithm  seems  to  be  much  better  than  the  previous  ones,  since  the  required 
time  to  solve  equations  from  the  AES  is  much  smaller  [5], 
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III.  COMPUTATIONAL  EXPERIMENTS 


A.  METHODOLOGY 

This  chapter  will  explain  how  the  MRHS  method  was  applied  by  using 
different  equations  that  were  generated  with  the  code  shown  in  Appendix  B.  The 
parameters  of  this  small-scale  variant  of  AES  were  determined:  the  number  of 
rounds  (0<Nr<10),  the  number  of  the  rows  and  the  columns  in  the  rectangular 
arrangement  of  the  input  (0<Nb<4),  and  the  number  of  bits  in  each  word  (2,  4,  8). 
In  addition,  a  random  plaintext  and  an  initial  key  were  chosen,  which  have 
lengths  that  correspond  to  the  number  of  columns  and  rows  of  the  states  and 
with  the  number  of  bits  in  the  words  as  well.  These  texts  are  expressed  in 
hexadecimal  notation.  After  that,  the  MRHS  Linear  Equations  method  was  used 
to  try  to  solve  this  system. 

Initially,  the  agreeing  procedure  was  implemented  to  see  if  some  linear 
equations  could  be  extracted  just  by  applying  this  procedure.  If  the  number  of 
linear  equations  that  resulted  was  equal  to  the  number  of  variables  in  the  system, 
then  the  procedure  would  be  stopped  since  there  would  be  a  unique  solution  for 
the  system.  If  the  number  of  linear  equations  extracted  from  the  agreeing 
procedure  was  less  than  the  number  of  variables,  the  gluing  procedure  was 
performed,  in  which  the  maximum  number  of  Right  Hand  Sides  (RHS)  in  the 
glued  equations  is  specified.  Then  another  attempt  was  made  to  extract  some 
linear  equations.  In  order  to  have  a  unique  solution,  it  is  necessary  to  have  a 
sum  of  linear  equations  equal  to  the  number  of  variables. 

B.  RESULTS 

Some  simple  examples  of  the  procedure  described  above  are  shown 

below. 
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1.  Example  1 


In  this  example,  a  system  was  created  with  two  rounds,  two  rows  and  two 
columns,  in  the  rectangular  arrangement  of  the  input  and  2-bit  words.  The 
system  has  also  a  plaintext  fa  (in  hexadecimal  notation),  which  corresponds  to 
11111010  in  binary  notation,  and  a  key  ea  (in  hexadecimal  notation),  which 
corresponds  to  11101010  in  binary  notation.  The  system  of  equations  appears 
below,  starting  with  the  number  of  bit  variables  and  the  number  of  MRHS 
equations.  Each  MRHS  equation  appears  as  the  number  of  rows,  the  number  of 
RHS,  the  rows  of  Aj  and  the  columns  of  Li.  The  first  four  symbols  are  from  the  key 
schedule,  the  next  four  from  round  1 ,  and  the  last  four  from  round  2. 

24  12 

4  4 

000000100000000000000000 
000000010000000000000000 
100000001000000000000000 
010000000100000000000000 
0011 
0110 
1000 
1101 
4  4 

000010000000000000000000 
000001000000000000000000 
001000000010000000000000 
000100000001000000000000 
0010 
0111 
1001 
1100 
4  4 

000000100010000000000000 

000000010001000000000000 

000000001000100000000000 


000000000100010000000000 
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0000 


000010001000000000000000 

000001000100000000000000 

000000000010001000000000 


000000000001000100000000 


100000000000000000000000 

010000000000000000000000 

000000000000000010000000 

000000000000000001000000 


000000100000000000000000 

000000010000000000000000 

000000000000000000100000 

000000000000000000010000 


000010000000000000000000 

000001000000000000000000 


000000000000000000001000 


000000000000000000000100 
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0010 


001000000000000000000000 

000100000000000000000000 

000000000000000000000010 


000000000000000000000001 


000000001000000001110000 

000000000100000011100000 

000000000000100000000000 

000000000000010000000000 


000000100010000000001101 

000000010001000000001011 

000000000000001000000000 

000000000000000100000000 


000010001000000000000111 

000001000100000000001110 


000010001000100000000000 


000001000100010000000000 
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0010 
0111 
1001 
1100 
4  4 

000000000010000011010000 

000000000001000010110000 

000000100010001000000000 

000000010001000100000000 

0000 

0101 

1011 

1110 

now:  neq  =  12  ;  nlink  =  0  ;  linbank  =  0 

Let  us  explain  the  above  shown  result.  The  “neq”  is  the  number  of  MRHS 
Linear  Equations  we  have  up  to  that  point.  The  “linbank”  (linear  bank)  is  the 
notation,  which  shows  the  number  of  the  ordinary  linear  equations  that  have 
extracted  after  the  application  of  the  agreeing  or  gluing  procedures. 

The  above  system  has  at  the  beginning  24  variables  combined  in  12 
MRHS  linear  equations.  Before  the  application  of  the  agreeing  procedure,  we 
cannot  extract  any  ordinary  linear  equation.  That  is  why  we  see  in  the  above 
result  that  we  have  0  number  of  linear  equation  in  our  linear  bank  (linbank). 

This  system  of  equations  was  used  as  input  in  the  code  “mrhs,”  which  is 
shown  in  Appendix  B.  After  the  agreeing  procedure,  we  conclude  with  zero 
MRHS  linear  equations  and  22  new  ordinary  linear  equations,  which  appear 
below. 


0 

+ 

8 

+ 

20 

+ 

21 

+ 

23  = 

0 

8 

+ 

9 

+ 

10 

+ 

21 

+ 

22  = 

0 

21  = 

0 

2 

+ 

10 

+ 

20  = 

0 

1 

+ 

2 

+ 

9 

+ 

10 

+ 

19  = 

0 

0 

+ 

8 

+ 

18  = 

0 
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0 

+ 

2 

+ 

8 

+ 

10 

+ 

17 

= 

0 

0 

+ 

8 

+ 

16 

= 

0 

15 

= 

0 

2 

+ 

8 

+ 

9 

+ 

10 

+ 

14 

= 

0 

0 

+ 

2 

+ 

8 

+ 

10 

+ 

13 

= 

0 

1 

+ 

2 

+ 

9 

+ 

10 

+ 

12 

= 

0 

1 

+ 

2 

+ 

9 

+ 

10 

+ 

11 

= 

0 

2 

+ 

8 

+ 

10 

= 

0 

0 

+ 

1 

+ 

2 

+ 

8 

+ 

9 

= 

0 

8 

= 

1 

2 

+ 

7 

= 

0 

0 

+ 

6 

= 

0 

5  =  0 
4  =  0 

0  +  2  +  3  =  0 

1  =  0 

now:  neq  =  0  ;  nlink  =  0  ;  linbank  =  22 

The  numbers  that  appear  in  the  above  system  of  ordinary  linear  equations 
represent  the  number  of  the  variable  in  our  initial  system,  i.e.,  0. 

The  matrix  form  of  the  above  ordinary  linear  equations  is  shown  below: 


001100000000000000000001 

001000000000000000000010 

000011000000000000000100 

000010000000000000001000 

000000110000000000010000 

000000100000000000100000 

110000000000000001000000 

100000000000000010000000 

000011001101000100000000 

000010001010001000000000 

000000110111010000000000 

000000101010100000000000 

000111000001000000000000 

001010000010000000000000 

110000011100000000000000 
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100000101000000000000000 

111010110000000000000000 

001101100000000000000000 

111111000000000000000000 

010010000000000000000000 

111100000000000000000000 

010000000000000000000000 

0101010101000101000000 

The  last  row  represents  the  Right  Hand  Side,  which  has  only  22  elements, 
meaning  that  there  are  two  free  variables.  The  same  thing  can  be  easily 
observed  from  the  reduced  matrix  of  the  coefficients  of  the  variables.  That 
particular  system  cannot  have  a  unique  solution.  The  possible  solutions  of  this 
system  appear  below. 

010000000000000000000000  0 
101100000000000000000000  0 
000010000000000000000000  0 
000001000000000000000000  0 
100000100000000000000000  0 
001000010000000000000000  0 
000000001000000000000000  1 
101000000100000000000000  1 
001000000010000000000000  1 
101000000001000000000000  0 
101000000000100000000000  0 
100000000000010000000000  0 
001000000000001000000000  1 
000000000000000100000000  0 
100000000000000010000000  1 
100000000000000001000000  0 
100000000000000000100000  1 
101000000000000000010000  0 
000000000000000000001000  1 
000000000000000000000100  0 
001000000000000000000010  1 
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100000000000000000000001  0 


free : 

0  2 


solutions : 

?0??00??1??????0????10?? 

000000001110001010101010 

100100101011111001011011 

001100011001100010111000 

101000111100010001001001 

Therefore,  the  variables  xo  and  X2  are  free  variables,  as  it  is  appeared  in 
the  code  above.  Additionally,  we  have  a  depiction  of  the  four  different  solutions  at 
the  end,  which  means  that  this  particular  system  has  four  different  keys  for  one 
plaintext  and  ciphertext.  In  that  example,  where  we  have  a  small-scale  variant  of 
the  AES  with  two  rounds,  two  rows  and  two  columns  in  the  rectangular 
arrangement  of  the  input  and  two-bit  words,  the  algorithm  cannot  give  a  unique 
solution. 

This  is  a  surprising  result,  as  we  mentioned  in  the  introduction.  We  have  a 
specific  plaintext,  which  is  encrypted  by  using  a  specific  key.  In  our  attempt  to 
break  the  algorithm  and  to  extract  the  initial  plaintext  from  the  ciphertext,  we 
expect  to  recover  the  initial  key.  Instead,  solving  the  system  of  the  MRHS  Linear 
Equations,  we  found  four  different  keys  that  can  decrypt  that  ciphertext  to  get  the 
correct  plaintext. 

One  thing  that  one  could  suppose  is  that  some  other  linear  equations 
could  be  extracted  in  order  to  solve  this  problem,  by  applying  a  different  plaintext 
with  the  same  key.  At  first  glance,  it  can  be  seen  that  the  ordinary  linear 
equations  are  different  from  the  previous  ones.  However,  if  the  matrix  of  the 
coefficients  of  the  variables  is  reduced,  we  conclude  exactly  in  the  same  set  of 
linear  equations  as  before.  Therefore,  such  a  system  has  a  non-unique  solution. 
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2.  Example  2 


The  second  example  is  one  with  a  unique  solution.  The  parameters  of  this 
system  are  one  round,  two  rows,  and  two  columns  in  the  rectangular 
arrangement  of  the  input  and  four  bit  words. 

The  system,  which  appears  below,  has  24  variables  and  six  MRHS  linear 
equations.  Therefore,  in  order  to  end  up  with  a  unique  solution,  it  is  necessary  to 
find  24  ordinary  linear  equations. 

24  6 

8  16 

000001000000000000000000 
000100000000000000000000 
001000000000000000000000 
000110000000000000000000 
010000000000000000000000 
011111010000000000000000 
011011100000000000000000 
111011100000000000000000 
10000000 
01000000 
10100000 
11110000 
00001000 
11001100 
10101100 
01101100 
00010110 
11000101 
10110101 
11100110 
00001101 
01011010 
00101101 
01101110 
8  16 
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000000000100000000000000 
000000010000000000000000 
000000100000000000000000 
000000011000000000000000 
000001000000000000000000 
000000110001000000000000 
000000100010000000000000 
000010000000000000000000 
10000000 
01000000 
10100000 
11110000 
00001000 
11001100 
10101010 
01101010 
00010111 
11000011 
10110011 
11100001 
00001101 
01011011 
00101101 
01101001 
8  16 

000000000000001000000000 

000000000001000000000000 

000000000010000000000000 

000000000001100000000000 

000000000100000000000000 

000000000011000100000000 

000000000100110000000000 

000000001000000000000000 

10000000 

11000000 

00100000 
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01110000 
10001000 
11001100 
10101010 
11101010 
00010101 
01000011 
10110011 
01100011 
10001101 
01011001 
00101101 
01101011 
8  16 

100000000000000000000000 

100100000000000000000000 

101000000000000000000000 

100100000000000010000000 

110000000000000000000000 

001100000000000000010000 

001000000000000000100000 

011000000000000011100000 

10000000 

11000000 

10100000 

11110000 

10001000 

11001100 

10101010 

11101010 

01101111 

00111010 

01001010 

00011001 

01110100 

00100011 
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01010100 


00010001 
8  16 

000000000000000000000100 
000000010000000000010000 
000000100000000000100000 
000000010000000000011000 
000001000000000001000000 
000000110000000000110001 
000001100000000001101110 
000011100000000011101110 
10000000 
01000000 
10100000 
11110000 
00001000 
11001100 
10101000 
01101000 
00010110 
11000001 
10110001 
11100010 
00001101 
01011010 
00101101 
01101010 
8  16 

000000000000010000000000 

000000000001000000000001 

000000000010000000000010 

000000000001100000000001 

000000000100000000000100 

000000000011000100000011 

000000000110111000000110 

000000001110111000001110 
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10000000 


01000000 

10100000 

11110000 

00001000 

11001100 

10101000 

01101000 

00010110 

11000001 

10110001 

11100010 

00001101 

01011010 

00101101 

01101010 

now:  neq  =  6  ;  nlink  =  0  ;  linbank  =  0 

In  this  case,  we  tried  to  extract  some  ordinary  linear  equations  directly 
from  the  initial  system  of  the  six  MRHS  Linear  Equations,  without  any  result.  After 
the  agreeing  algorithm  is  applied,  there  is  still  no  linear  equation  in  the  linear 
bank.  As  we  said  before,  “linbank”  (Linear  bank)  is  the  name  of  the  place  in  the 
computer  program  where  all  the  ordinary  linear  equations,  that  are  extracted  from 
the  different  procedures,  are  saved.  This  means  that  many  times  the  agreeing 
procedure  has  no  result,  if  it  is  applied  alone  in  such  a  system  of  MRHS  Linear 
Equations.  Here  is  where  the  gluing  procedure  is  performed.  In  this  particular 
example,  we  observe  that  exactly  24  different  ordinary  linear  equations  are 
extracted,  which  are  shown  below. 


4  +  23  =  0 

4  +  22  =  0 

4  +  21  =  0 

4  +  20  =  0 
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19  =  0 
17  +  18  =  0 

17  =  0 
16  =  0 
4  +  15  =  0 

14  =  0 
4  +  13  =  0 

12  =  0 
4  +  11  =  0 

10  =  0 
4  +  9  =  0 

4  +  8  =  0 

4  +  6  +  7  =  0 

4  +  6  =  0 

5  =  0 

3  +  4  =  0 

3  =  1 

2  =  0 

1  =  0 

0  =  0 

now:  neq  =  0  ;  nlink  =  0  ;  linbank  =  24 

As  we  mentioned  before,  we  conclude  in  zero  MRHS  Linear  Equations 
and  24  ordinary  linear  equations  (“linbank=24”). 

In  a  matrix  form,  these  are: 

24  0 

linear  eqs  (24 )  : 

24  1 

000010000000000000000001 

000010000000000000000010 

000010000000000000000100 

000010000000000000001000 

000000000000000000010000 

000000000000000001100000 

000000000000000001000000 

000000000000000010000000 
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000010000000000100000000 


000000000000001000000000 

000010000000010000000000 

000000000000100000000000 

000010000001000000000000 

000000000010000000000000 

000010000100000000000000 

000010001000000000000000 

000010110000000000000000 

000010100000000000000000 

000001000000000000000000 

000110000000000000000000 

000100000000000000000000 

001000000000000000000000 

010000000000000000000000 

100000000000000000000000 

000000000000000000001000 

In  addition,  the  unique  solution  of  this  system  is: 

100000000000000000000000  0 
010000000000000000000000  0 
001000000000000000000000  0 
000100000000000000000000  1 
000010000000000000000000  1 
000001000000000000000000  0 
000000100000000000000000  1 
000000010000000000000000  0 
000000001000000000000000  1 
000000000100000000000000  1 
000000000010000000000000  0 
000000000001000000000000  1 
000000000000100000000000  0 
000000000000010000000000  1 
000000000000001000000000  0 
000000000000000100000000  1 
000000000000000010000000  0 
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000000000000000001000000  0 


000000000000000000100000  0 
000000000000000000010000  0 
000000000000000000001000  1 
000000000000000000000100  1 
000000000000000000000010  1 
000000000000000000000001  1 


free : 


solutions : 

000110101101010100001111 

There  are  no  free  variables.  Therefore,  it  is  a  unique  solution,  as  it 
appears  above. 

From  the  two  examples  above,  we  saw  how  the  algebraic  attack  method 
of  MRHS  Linear  Equations  works.  The  many  experiments  that  were  tested  for 
this  thesis  were  mostly  successful  for  a  small  variant  of  the  AES.  The  parameters 
of  the  system  that  were  used  were  at  most  four  rounds,  with  four  rows  and 
columns  and  8-bit  words  instead  of  10  rounds  that  the  original  AES  algorithm 
uses. 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  work  has  shown  that  the  new  method  of  breaking  the  Rijndael 
algorithm  with  the  algebraic  representation  of  a  system  with  Multiple  Right  Hand 
Sides  linear  equations  is  quite  effective.  It  is  comparatively  effective  because  by 
using  other  algebraic  methods,  such  as  the  Buchberger’s  and  the  F4  -  F5 
algorithms,  it  was  impossible  to  find  a  solution  for  some  very  simple  systems.  In 
particular,  a  small  version  of  the  AES  with  5  rounds  and  1  row  and  1  column  in 
the  rectangular  arrangement  of  the  input,  could  not  be  solved  with  the  older 
algebraic  methods  [5],  Moreover,  a  system  with  4  rounds  and  1  row  and  1 
column  in  the  rectangular  arrangement  of  the  input,  took  20286.18  seconds  to  be 
solved  by  the  standard  algebraic  methods  [5],  The  MRHS  method  solved  this 
problem  in  only  0.032  seconds  [5], 

In  Chapter  II,  we  saw  that  the  new  method  seems  to  be  much  better  than 
the  other  algebraic  methods  in  terms  of  complexity.  As  mentioned  before,  due  to 
the  different  parameters  that  this  algorithm  can  have,  (different  number  of 
rounds,  rows  and  columns)  it  is  very  complicated  to  specify  an  exact  complexity. 
According  to  the  real  parameters  of  the  AES,  our  test  cannot  be  considered  very 
realistic,  since  the  values  of  the  parameters  of  the  original  AES  are  bigger. 
Nevertheless,  by  comparing  the  complexities  in  terms  of  the  time  consumed  for  a 
small-scale  variant  to  be  solved,  we  conclude  that  it  is  worthwhile  to  further 
develop  and  improve  this  method. 

The  AES,  as  we  explained  in  the  beginning,  is  a  standard  in  order  to 
secure  federal  information.  With  MRHS  Linear  Equations,  we  concluded  that  we 
can  break  a  small-scale  variant  of  this  system.  We  cannot  say  that  the  AES  is  at 
risk,  but  we  can  say  for  sure  that  this  new  method  gives  new  perspectives  in  the 
field  of  cryptanalysis  that  may  put  the  AES  algorithm  at  risk. 

For  example,  one  new  detail  in  this  thesis,  is  the  application  of  our  codes 
in  many  different  experiments  with  different  parameters.  In  particular,  instead  of 
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the  8-bit  field  that  professors  Havard  Raddum  and  Igor  Semaev  used,  we 
executed  some  experiments  in  the  2-bit  and  4-bit  fields,  that  is  original.  From 
that,  we  discovered  a  result  very  important  in  terms  of  cryptanalysis.  There  are 
some  cases  where  (a  small  variant  of)  the  Rijndael  algorithm  can  have  multiple 
solutions.  We  demonstrated  that  in  the  2-bit  field.  However,  it  may  be  valid 
though  more  difficult  of  course,  even  in  the  8-bit  field,  which  might  reduce  the 
effectiveness  of  the  algorithm.  This  means  that  we  may  have  to  be  cautious 
when  we  choose  the  encryption  keys,  because  multiple  decryption  keys  would 
make  it  easier  to  break  the  algorithm  and  further  to  decrypt  and  reveal  important 
federal  information. 

In  conclusion,  our  first  goal  was  to  create  the  code,  which  could  solve  a 
small-scale  variant  of  the  AES.  Based  on  the  codes  of  Professors  Raddum  and 
Semaev,  we  reached  this  point.  In  a  few  cases,  a  unique  solution  of  a  small-scale 
variant  of  the  AES  resulted.  However,  the  small  AES  algorithm  has  some 
weaknesses,  which  consist  of  the  fact  that  some  systems  cannot  uniquely  be 
solved,  even  if  different  plaintexts  with  the  same  initial  key  are  applied. 

The  complexity  of  this  algorithm  is  the  area  that  needs  to  be  improved  in 
future  research.  The  reduction  of  the  complexity  by  a  significant  factor  could 
upgrade  the  MRHS  algorithm  to  a  very  effective  algebraic  attack  method  against 
new  very  strong  cryptographic  algorithms.  This  reduction  could  be  achieved  by 
reducing  either  the  number  of  the  agreeings  between  the  symbols  or  the 
numbers  of  gluings  which  adds  more  complexity  in  the  whole  algorithm.  One  way 
to  achieve  that  could  be  the  simultaneous  agreeing  of  more  than  two  symbols, 
which  may  require  a  smaller  number  of  necessary  gluings  for  the  extraction  of 
the  ordinary  linear  equations.  In  that  way,  we  could  break  an  AES  algorithm  with 
parameters  of  greater  value. 
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APPENDIX  A 


(These  codes  are  only  intended  for  experiments  with  the  MRHS  method, 
and  neither  their  author  nor  the  author  of  this  thesis  takes  any  responsibility  for 
their  use) 

1.  BASIS. h  (by  Havard  Raddum  and  Igor  Semaev) 

#include  <stdio.h> 

#include  <stdlib.h> 

#include  <math.h> 


typedef  unsigned  int  u32; 
typedef  unsigned  char  u8; 

#def ine  MAXINT  2147483647 
#def ine  EPSILON  0.001 

#def ine  log2 (x)  (  1.442695040888963407359924681  *  log(x)  )  //  ! ! 

changed  for  bcc  -  DC 

int  NVAR,  NWORDS; 


u8 

weight  [256]={0, 1,1, 2, 1,2, 2, 3, 1,2, 2, 3, 2, 3, 3, 4, 1,2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 

4, 4, 5, 1,2, 2, 3, 2, 3, 3, 4, 

2. 3. 3. 4. 3. 4. 4. 5. 2. 3. 3. 4. 3. 4. 4. 5. 3. 4. 4. 5. 4. 5. 5. 6. 1.2. 2. 3. 2. 3. 3. 4. 2 
3343445 

,  J,  h,  J,  h,  h, 

2. 3. 3. 4. 3. 4. 4. 5. 3. 4. 4. 5. 4. 5. 5. 6. 2. 3. 3. 4. 3. 4. 4. 5. 3. 4. 4. 5. 4. 5. 5. 6. 3 

,4, 4, 5, 4, 5, 5, 6, 

4. 5. 5. 6. 5. 6. 6. 7. 1.2. 2. 3. 2. 3. 3. 4. 2. 3. 3. 4. 3. 4. 4. 5. 2. 3. 3. 4. 3. 4. 4. 5. 3 

,4, 4, 5, 4, 5, 5, 6, 

2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 2 

3343445 


3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 3, 4, 4, 5, 4, 5, 5, 6, 4 

,5, 5, 6, 5, 6, 6, 7, 

4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7,7,8}; 


struct  bitVector{ 
u32  *v; 

int  length,  wl; 

}  ; 


struct  eqSymbol{ 

int  nlin,  nrhs,  nw,  nn,  eqnr,  originalNRHS ; 
u32  * *A,  **b,  **coverID; 
u8  *RHSexists,  delSinceExtract; 
struct  linkSymbol  **link; 
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In¬ 


struct  linkSymbol{ 

int  ncells,  *nCover[2],  linknr; 
struct  eqSymbol  *nlist[2]; 
struct  bitVector  *cellID; 
u8  *cellExists; 

In¬ 


struct  system{ 
int  neq,  nlink; 
struct  eqSymbol  *E,  *linbank; 
struct  linkSymbol  *L; 

}  ; 

int  ww(u32  *word,  int  n) { 
int  i,  j,  w=0; 

for(i=0;  i<n;  ++i){ 

w+=weight [word [ i ] &0xf f ] ; 
w+=weight [ (word [ i ] >>  8 ) &0xf f ] ; 
w+=weight [ (word [i] >>16) &0xff]  ; 
w+=weight [word[i] >>24] ; 

} 

return  w; 


double  averageNRHS ( struct  system  *S) { 
int  i,  totnrhs=0; 


if ( S->neq==0 ) 
return  0.0; 

for(i=0;  i<S->neq;  ++i) 
totnrhs+=S->E [ i ] . nrhs; 

return  (double) (totnrhs) / (double) (S->neq) 


int  smallestSetBit (u32  *M,  int  new) { 
int  j=0,  m=0; 
u32  w; 

while(j<ncw  &&  !M[j]){ 

j++; 
m+=32 ; 

} 

if ( j  ==ncw) 
return  -1; 
w=M [ j ] ; 

for ( j=4;  j  >=0 ;  --j )  { 

if  (!  (w&  (  (1«  (l«j  )  )  -1)  )  )  { 
m+=  (l«j  )  ; 
w >>=  ( 1«  j  )  ; 

} 
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} 


return  m; 


int  largestSetBit (u32  *M,  int  new) { 
int  j ,  m; 
u32  w; 

j=ncw-l; 
m= (ncw<<5 ) -1 ; 
while(j>=0  &&  !M[j]){ 

j--; 

m-=32 ; 

} 

if ( j= — 1) 

return  -1; 
w=M [ j ] ; 

for ( j=4;  j  >=0 ;  --j )  { 

if(!(w&(((l<<(l<<j))-l)AOxffffffff)))//no  set  bit  in  upper  remaning 

half 

m-=  (l«j  )  ; 
else 

w >>=  ( 1«  j  )  ; 

} 

return  m; 

} 

u8  vOGreaterThanVl ( struct  bitVector  vO,  struct  bitVector  vl){ 
int  i; 

if (vO . length ! =vl . length) printf ( " (vOGreaterThanVl ) Uncomparable  bit- 
strings,  vO  . length=%d,  vl . length=%d ! \n"  , 

vO  . length, vl . length) ; 

i=vO .wl-1; 

while (i>=0  &&  !(v0.v[i]Avl.v[i])) 

i--; 

if(i>=0  &&  vO . v [ i ] >vl . v [ i ] ) 
return  1; 
else 

return  0; 


u8  vOEqualVl ( struct  bitVector  vO,  struct  bitVector  vl){ 
int  i; 

if (vO . length ! =vl . length) printf ( " (vOEqualVl ) Uncomparable  bit- strings , 
vO . length=%d,  vl . length=%d ! \n" , 

vO . length, vl . length) ; 

i=vO . wl-1; 

while (i>=0  &&  !(v0.v[i]Avl.v[i])) 

i— ; 
if (i<0) 
return  1; 
else 

return  0; 
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} 


void  mergeSortBitVectors ( struct  bitVector  *vl,  int  nil){ 
int  i,  j,  t,  na,  nb; 
struct  bitVector  tmp,  *ml,  *vlb; 


if (nil==l ) 
return; 
if (nil==2 ) { 

if (vOGreaterThanVl (vl [0]  ,vl [1] ) )  { 
tmp=vl [ 0 ] ; 
vl [ 0 ] =vl [ 1 ] ; 
vl [ 1 ] =tmp; 

} 

return; 

} 

na=nil/ 2 ; 

nb=nil-na; 

vlb=vl+na; 

mergeSortBitVectors (vl,na) ; 
mergeSortBitVectors (vlb, nb) ; 

ml= (struct  bitVector  *) malloc (nil* sizeof ( struct  bitVector)); 
t=i= j  =0 ; 

while (i<na  &&  j<nb) { 

if (vOGreaterThanVl (vl [ i ] , vlb [ j ] ) ) 
ml [ t++ ] =vlb [ j  ++ ] ; 
else 

ml [t++] — v 1 [i++] ; 

} //lists  merged 
if ( i<na) { 

for ( j  =i ;  j  <na;  ++ j ) 
ml [t++] =vl [ j ] ; 

} 

else  { 

for(i=j;  i<nb;  ++i) 
ml [t++] =vlb [i] ; 

} //remainder  of  unfinished  list  copied 
for(i=0;  i<nil;  ++i) 
vl [ i ] =ml [ i ] ; 

//copied  back  to  vl 


void  printLinComb (u32  *1,  int  nvar) { 
int  i; 


for(i=0;  i<nvar;  ++i){ 

if  (1  [i»5]  &  (1«  (i&Oxlf )  )  ) 
printf("%3d  +  ",i); 

} 

printf ( "\b\b\n" ) ; 


void  printEquation ( struct  eqSymbol  *eq) { 

/*  Prints  the  equation  to  the  screen,  includes  up  to  16  righ-hand 
sides  */ 
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int  i,  j,  maxnv=0,  nv,  nr,  nwritten; 

if (eq->nrhs>16) 
nr=16; 

else 

nr=eq->nrhs ; 

for(i=0;  i<eq->nlin;  ++i){ 
nv=ww (eq->A [ i ] , NWORDS )  ; 
if (nv>maxnv) 
maxnv=nv; 

} 

printf ( "=========  Equation  %d  ================================\n" , eq- 

>eqnr) ; 

for(i=0;  i<eq->nlin;  ++i){ 
nv=ww (eq->A [ i ] , NWORDS ) ; 
for(j=0;  j<maxnv-nv;  ++ j ) 
printf ("  "); 

for(j=0;  j<NVAR;  ++ j ) { 

if  (eq->A[i]  [  j»5]  &  (1<<  (j&Oxlf)  )  ) 
printf ("  %3d  +",j); 

} 

printf ("\b=") ; 
nwritten=0 ; 

for(j=0;  j <eq->originalNRHS ;  ++ j ) { 
if (eq->RHSexists  [  j ] )  { 
if  (eq->b  [j]  [i»5]&(l«(i&0xlf))) 
printf ("  1  I"); 
else 

printf ("  0  !"); 
nwritten++; 
if (nwritten==nr ) 

j  =eq->originalNRHS  ; 

} 

} 

printf ("\b\n\n") ; 

} 

printf ("%d  right  hand  sides  in  total\n", eq->nrhs) ; 

printf (" \nLinked  to  %d  other  equations\n" , eq->nn) ; 

printf ("RHS  index  -  "); 

nwritten=0 ; 

for(i=0;  i<eq->originalNRHS ;  ++i){ 
if (eq->RHSexists [i] ) { 
printf ("  %2d  ", i) ; 
nwritten++; 
if (nwritten==nr ) 
i=eq->originalNRHS  ; 

} 

} 

printf ( " \ncoverID" ) ; 

for(i=0;  i<eq->nn;  ++i){ 

printf ("\nLink  %2d  ",i); 

nwritten=0  ; 

for(j=0;  j <eq->originalNRHS ;  ++ j  )  { 
if (eq->RHSexists [ j ] ) { 
printf ("  %2d  " , eq->coverID [ i ]  [ j ]  )  ; 
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nwritten++; 
if (nwritten==nr ) 

j  =eq->originalNRHS  ; 

} 

} 

} 

printf ( " \n=============== 


void  printLink ( struct  linkSymbol  *ls){ 

/*  Prints  the  link  to  the  screen  */ 
int  i,  cO,  cl; 

printf ("==============  Link  %d  ========================\n" , ls- 

>linknr)  ; 

printf ("Link  with  %d  cells\n" , ls->ncells) ; 
if ( ls->ncells<=32 ) { 
printf ("\n  "); 

for(i=0;  i<ls->ncells ;  ++i) 
printf ("  %2d  ",  i)  ; 
printf (" \nExists  "); 
for(i=0;  i<ls->ncells ;  ++i) { 
if ( ls->cellExists [i]  ) 
printf ("  x  "); 
else 

printf  ("  "); 

} 

printf ( "\nnCoverO" ) ; 
for(i=0;  i<ls->ncells ;  ++i) 

printf ("  %2d  " , ls->nCover [0] [i] ) ; 
printf ( " \nnCover 1 " ) ; 
for(i=0;  i<ls->ncells ;  ++i) 

printf ("  %2d  " , ls->nCover [ 1 ]  [ i ]  )  ; 

} 

else  { 

cO=ls->nCover  [ 0 ]  [0]; 
cl=ls->nCover [ 1 ]  [0]; 
for(i=l;  i<ls->ncells ;  ++i){ 

if (cO ! =ls->nCover [ 0 ] [ i ]  | |  cl ! =ls->nCover [ 1 ] [ i ] ) { 

printf ( " \nUnbalanced  link\n") ; 
i=ls->ncells ; 

} 

} 

} 

printf (" \n\nLinks  together  equations  %d  and  %d  of  dimensions  %d 
%d\n" , 

ls->nlist [0] ->eqnr, ls->nlist [1] ->eqnr, ls->nlist [0] ->nlin, Is 
>nlist [1] ->nlin) ; 

printf ( "===================================================\n" ) 


void  printLinEquation ( struct  eqSymbol  *eq,  int  start,  int  stopp) { 
int  i,  j,  nv,  maxnv=0;  //  ! !  added  init  ,  maxnv; 

for(i=start;  ikstopp;  ++i) { 


\n" ) 


and 
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nv=ww (eq->A [ i ] , NWORDS )  ; 
if (nv>maxnv) 
maxnv=nv; 

} 

printf ( "=========  Linear  Equation  =====================\n" , eq->eqnr) 

for(i=start;  i<stopp;  ++i) { 
nv=ww (eq->A [ i ] , NWORDS )  ; 
for(j=0;  j<maxnv-nv;  ++ j ) 
printf ("  "); 

for(j=0;  j<NVAR;  ++ j ) { 

if (eq->A[i] [j>>5] & (1<< (j&Oxlf) ) ) 
printf ("  %3d  +",j); 

} 

printf ("\b=") ; 

if (eq->b [0] [i>>5] & (1<< (i&Oxlf) ) ) 
printf ("  l\n") ; 
else 

printf ("  0\n"); 

} 

} 

void  deleteEquation ( struct  eqSymbol  *eq) { 

/*  Frees  all  allocated  memory  in  eq  not  associated  to  links.  */ 
int  i; 


for(i=0;  i<eq->originalNRHS ;  ++i){ 
if (eq->RHSexists [i] ) 
free (eq->b [ i ] ) ; 

} 

free (eq->b) ; 

for(i=0;  i<eq->nlin;  ++i) 
free (eq->A [ i ] ) ; 
free (eq->A) ; 
free (eq->RHSexists ) ; 


void  deleteSystem ( struct  system  *S) { 
int  i; 


for(i=0;  i<S->neq;  ++i) 
deleteEquation ( S->E+i ) ; 
deleteEquation ( S->linbank)  ; 


void  Uxb(u32  **U,  int  nr,  int  nc,  u32  *b,  u32  *x) { 

/*  computes  U  times  b,  where  U  has  nr  rows  and  nc  columns.  Stores 
result  in  x  */ 
int  i,  j,  new,  nrw; 
u32  *w; 

ncw= (nc+31) >>5; 
nrw= (nr+31) >>5; 

w— (u32  * ) malloc (new* sizeof (u32 ) ) ; 
for(i=0;  i<nrw;  ++i) 
x [ i ] =0 ; 


55 


for(i=0;  i<nr;  ++i){ 
for(j=0;  j<ncw;  ++ j ) 
w  [  j ] =U [ i ]  [ j ] &b  [  j  ]  ; 
if (ww (w, new) &1 ) 

x[i>>5]  |  =  (l«(i&Oxlf)  )  ; 

} 

free (w)  ; 


void  UAX(u32  **U,  u32  **A,  u32  **nyA,  int  nl,  int  new) { 

/*  Multiplies  the  matrices  U  and  A,  stores  result  in  nyA  */ 
int  i,  j,  k,  *bpos,  hw; 

bpos= (int  * ) malloc (nl*sizeof (int) ) ; 
for(i=0;  i<nl;  ++i){ 
hw=  0  ; 

for ( j=0;  j  <nl ;  ++j )  { 

if  (U  [i]  [  j»5]  &  (1«  (j&Oxlf)  )  ) 
bpos [hw++] =j ; 

} //found  positions  where  U[i]  has  1-bits 
for(j=0;  j<ncw;  ++ j ) { 
for(k=0;  k<hw;  ++k) 
nyA [ i ] [ j ] A=A [bpos [k] ] [ j ] ; 

} 

} 

} 

int  rank(u32  **M,  int  nr,  int  nc) { 

/*  Returns  the  rank  of  M  */ 

int  rank,  i,  j,  k,  sb,  minsb,  r,  new; 

u32  **cM,  *tmp; 

ncw= (nc+31) >>5; 

cM= (u32  **) malloc (nr*sizeof (u32) ) ; 
for(i=0;  i<nr;  ++i){ 

cM [ i ]  =  (u32  * ) malloc (new* sizeof (u32 ) )  ; 
for(j=0;  j<ncw;  ++ j  ) 
cM[i]  [  j  ]  =M[i]  [j  ]  ; 

} 

for(i=0;  i<nr;  ++i) { 
minsb=nc; 

for ( j  =i ;  j  <nr ;  ++ j )  { 

sb=smallestSetBit (cM[ j ] , new) ; 
if(sb!=-l  &&  sb<minsb) { 
r=j  ; 

minsb=sb; 

} 

} 

if (minsb==nc) { 

for ( j  =0 ;  j  <nr ;  ++ j ) 
free (cM[ j ] ) ; 
free (cM) ; 
return  i; 

} 

if (r>i) { 
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tmp=cM [ i ] ; 
cM[i] =cM[r] ; 
cM [r] =tmp; 


} 

for  ( j=i+l;  j <nr ;  ++ j )  { 

if (cM[ j ] [minsb>>5] & (1<< (minsb&Oxlf ) ) ) { 
for(k=0;  k<ncw;  ++k) 
cM [ j ] [k] A=cM[i] [k] ; 

} 

} 

} 

for(i=0;  i<nr;  ++i) 
free (cM[i] ) ; 
free (cM) ; 
return  nr; 


int  computeU ( struct  eqSymbol  *eqO,  struct  eqSymbol  *eql,  u32  * 
/*  Computes  U  such  that  UM  is  triangularized,  where  M  is  the 
concatenation  of  the  matrices  in  eqO  and  eql .  Returns  the 
dimension  of  the  shared  subspaces  of  eqO  and  eql.  */ 
int  i,  j,  k,  b,  minb,  r,  enord,  uw,  rank,  sumnl; 
u32  *tmp,  enmask,  **M; 

sumnl=eqO->nlin+eql->nlin ; 

M=(u32  **) malloc ( sumnl*sizeof (u32  *)); 
for(i=0;  iksumnl;  ++i){ 

M[i] = (u32  *) malloc (NWORDS*sizeof (u32) ) ; 
if ( i<eqO->nlin) { 

for ( j  =0 ;  j  <NWORDS ;  ++ j  ) 

M[i]  [  j ] =eqO->A [ i ]  [ j ]  ; 

} 

else  { 

for ( j  =0 ;  j  <NWORDS ;  ++ j  ) 

M[i] [ j ] =eql->A[i-eqO->nlin] [j]; 

} 

} 

rank=sumnl ; 

uw= (sumnl+31) >>5; 

for(i=0;  i<sumnl;  ++i){ 

U[i]=(u32  *)calloc(uw,sizeof (u32) ) ; 

U[i]  [i>>5]  |  =  (1«  (i&Oxlf )  )  ; 

} 

//U  is  identity-matrix 
for(i=0;  iksumnl;  ++i) { 
minb=MAXINT; 
for(j=i;  jksumnl;  ++ j ) { 

b=smallestSetBit (M [ j ] , NWORDS ) ; 
if(b!=-l  &&  b<minb) { 
r=j  ; 
minb=b; 

if (minb==i) //no  need  to  search  further 
j  =sumnl ; 

} 

} 


U)  { 
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if (minb==MAXINT) { //only  all-zero  rows  remaining 
rank=i ; 
i=sumnl ; 

} 

else  { 

if (r>i) { //need  to  swap  rows 
tmp=M [ i ] ; 

M[i] =M[r] ; 

M [r] =tmp; 
tmp=U [ i ] ; 

U[i]=U[r]  ; 

U [r] =tmp; 

} 

enmask=l«  (minb&Oxlf )  ; 
enord=minb»5  ; 

for(j=i+l;  j<sumnl;  ++ j ){ //making  0's  under  leading  1 
if (M[ j ] [enord] Senmask) { 
for ( k=0 ;  k<NWORDS;  ++k) 

M  [  j ]  [k] A=M[i]  [k]  ; 
for(k=0;  k<uw;  ++k) 

U [ j ] [k] A=U[i] [k] ; 

} 

} 

} 

} 

for(i=0;  iksumnl;  ++i) 
free (M [ i ] ) ; 
free (M) ; 

return  sumnl-rank; 


u8  checkSolution ( struct  system  *S,  u8  *val) { 

/*  Returns  1  if  val  is  a  solution  to  S,  0  if  not. 

Prints  out  the  equation  number  of  equations  not  satisfied.  val[i] 
is  assigned  to  variable  i  .*/ 
int  i,  j; 
u32  *propRHS; 
struct  eqSymbol  *eq; 
u8  bit,  funnet,  rv=l; 


propRHS= (u32  *)calloc (NWORDS , sizeof (u32) ) ; 
for(eq=S->E;  eq<S->E+S->neq;  ++eq) { 
for(i=0;  i<eq->nlin;  ++i){ 
bit=0 ; 

for(j=0;  j<NVAR;  ++ j ) { 

if  (eq->A[i]  [  j»5]  &  (1«  ( j&Oxlf )  )  &&  val[j]) 

bit A=1 ; 

} 

if (bit) 

propRHS [i>>5] | =1<< (i&Oxlf ) ; 

}//propRHS  maa  vaere  blant  eq->b 
for(i=0;  i<eq->originalNRHS ;  ++i){ 
if (eq->RHSexists [i] ) { 
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funnet=l ; 

for(j=0;  j<eq->nw;  ++ j ) { 

if (propRHS [ j ] Aeq->b  [ i ]  [  j ] ) 
funnet=0 ; 

} 

if (funnet) 

i=eq->originalNRHS+l ; 

} 

} 

if ( i==eq->originalNRHS ) { //eq  ikke  tilf redsstilt 

printf ( "Equation  %d  not  satisfied,  nlin=%d,  nrhs=%d\n" , eq- 
>eqnr , eq->nlin, eq->nrhs ) ; 

/ /printEquation (eq)  ; 
rv=  0  ; 

} 

for(i=0;  i<eq->nw;  ++i) 
propRHS [ i ] =0 ; 

} 

return  rv; 

} 

u8  deleteRHS ( struct  eqSymbol  *eq,  int  RHSindex) { 
int  i,  minindex,  ci; 
struct  linkSymbol  *1S; 

if (eq->RHSexists [RHSindex] ==1 ) { 
eq->RHSexists [RHSindex] =0 ; 
free (eq->b [RHSindex] ) ; 
eq->nrhs-~ ; 
if (eq->nrhs==0 ) 
return  0; 

for(i=0;  i<eq->nn;  ++i){ 
lS=eq->link [ i ] ; 
ci=eq->coverID [ i ] [RHSindex]; 
if (lS->nlist [0] ==eq) 
minindex=0 ; 
else 

minindex=l ; 

lS->nCover [minindex] [ci]--; 

if (lS->nCover [minindex] [ci] <0) {print Link (IS) ;exit (0) ; } 

} 

eq->delSinceExtract=l ; 

} 

return  1; 


u8  trimEquation ( struct  eqSymbol  *eq) { 

/*  Makes  sure  that  A-matrix  in  eq  has  full  rank,  and  removes 

RHS  if  possible.  Returns  0  if  eq  can  not  be  satisfied,  returns  1 
if  A-matrix  had  full  rank,  returns  2  if  A-matrix  has  decreased.  */ 
u32  **U,  *x,  *tmp,  enmask; 
int  i,  j,  k,  b,  minb,  r,  enord,  uw,  rank; 
u8  rv;  //  ! !  removed  var:  ,  delocc=0; 

rank=eq->nlin; 
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uw= (eq->nlin+31 ) >>5; 

U=(u32  **) malloc (eq->nlin*sizeof (u32  *)); 
for(i=0;  i<eq->nlin;  ++i){ 

U[i]=(u32  *)calloc(uw,sizeof(u32) ) ; 

U[i]  [i»5]  |  =  (1«  (i&Oxlf )  )  ; 

} 

//U  is  identity-matrix 
for(i=0;  i<eq->nlin;  ++i){ 
minb=NVAR; 

for(j=i;  j<eq->nlin;  ++ j ) { 

b=smallestSetBit (eq->A [ j ]  ,  NWORDS )  ; 
if (b ! =— 1  &&  b<minb) { 
r=j  ; 
minb=b; 

} 

} 

if (minb==NVAR) { //only  all-zero  rows  remaining 
rank=i ; 
i=eq->nlin; 

} 

else  { 

if (r>i) { //need  to  swap  rows 
tmp=eq->A[i] ; 
eq->A[i] =eq->A[r] ; 
eq->A [r] =tmp; 
tmp=U [ i ] ; 

U[i]=U[r]  ; 

U [r] =tmp; 

} 

enmask=l<< (minb&Oxlf ) ; 
enord=minb>>5 ; 

for(j=i+l;  j<eq->nlin;  ++ j ){ //making  0's  under  leading  1 
if (eq->A [ j ] [enord] Senmask) { 
for ( k=0 ;  k<NWORDS;  ++k) 

eq->A[ j ] [k] A=eq->A[i] [k] ; 
for(k=0;  k<uw;  ++k) 

U  [  j ]  [k] A=U[i]  [k]  ; 

} 

} 

} 

}//U  is  computed 

x— (u32  * ) malloc (eq->nw*sizeof (u32  )  )  ; 
for(i=0;  i<eq->originalNRHS ;  ++i){ 
if (eq->RHSexists [i] ) { 

Uxb (U, eq->nlin, eq->nlin, eq->b [ i ] , x) ; //x  er  ny  RHS-vektor 

for(j=0;  j<eq->nw;  ++ j ) 

eq->b [ i ] [ j ] =x[ j ] ; 

for(j=rank;  j<eq->nlin;  ++ j ) { 

if (eq->b [i]  [ j>>5] & (1« ( j &0xlf )  )  )  { //inneholder  1-bit  i  0-rad 
omraade 

rv=deleteRHS (eq,  i)  ; 

//  ! !  removed  var:  delocc=l; 
if (rv==0) 
return  0; 
j=eq->nlin; 
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} 

} 

} 

}//oppdatert  alle  RHS,  og  sikret  at  de  er  gyldige 
if (eq->nlin==rank) 
return  1; 
else  { 

eq->nlin=rank; 

eq->nw= (eq->nlin+31 )  >>5; 

return  2; 

} 


void  depositLinComb ( struct  eqSymbol  *eq,  u32  *lc,  u8  rhv) { 

/*  Adds  lc=rhv  to  eq  (the  linear  bank) .  The  bank  must  be 

antitriangular ,  and  will  remain  so  after  addition  of  lc=rhv. 
lc  must  not  be  in  the  span  of  eq->A.  */ 
int  i,  j,  vnr,  insertrow=-l ,  libit; 
u32  *clc; 

clc= (u32  * ) malloc (NWORDS* si zeof (u32 ) ) ; 
for ( i=0 ;  i<NWORDS;  ++i) 
clc [i] =lc [i] ; 

//clc  er  kopi  og  kan  herjes  med 
if (eq->nlin==0 ) 
insertrow=0 ; 
i=0  ; 

while (insertrow==-l) { 

vnr=largestSetBit (clc,  NWORDS)  ; 

while ( i<eq->nlin  &&  vnrklargestSetBit (eq->A [ i ], NWORDS ) ) 
i++; 

if (i<eq->nlin) { 

libit=largestSetBit (eq->A [ i ] , NWORDS ) ; 
if (vnr>libit) 
insertrow=i ; 
else{ / /vnr==libit 
for ( j  =0 ;  j  <NWORDS ;  ++ j ) 
clc [ j ] A=eq->A[i] [ j ] ; 
if (eq->b[0] [i>>5] & (1<< (i&Oxlf) ) ) 
rhvA=l ; 

i++; 

} 

} 

else/ / i==eq->nlin 

insertrow=eq->nlin; 

} 

for (i=eq->nlin;  i>insertrow;  — i )  { 
eq->A[i] =eq->A[i-l] ; 

if  (  (  (eq->b  [0]  [i>>5]>>(i&0xlf))A(eq->b[0]  [  (i-1)  »5]  »  (  (i- 
1) &0xlf ) ) ) &1) 

eq->b [0] [i>>5]A=(l<<(i&0xlf)); 

}//ryddet  plass  til  ny  linear  ligning 
eq->A [ insertrow] =clc; 

if ( (u8) ( (eq->b [0] [ insert row>>5 ] >> (insertrow&Oxlf ) ) &1) Arhv) //feil  bit 
i  RHS 
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eq->b [0] [ insertrow>>5 ] A= (1« (insertrow&Oxlf ) ) ; 
eq->nlin++; 

} 

u8  substituteLinComb ( struct  system  *S,  u32  *lc,  u32  rhv) { 

/*  Largest  variable  in  lc  is  eliminated  using  the  equation  lc  =  rhv. 
Returns  0  if  system  becomes  inconsistent,  1  if  no  A-matrices  were 
reduced,  and  2  if  some  A-matrices  did  not  have  full  rank  after 
substitution.  */ 

int  i,  j,  enord,  maxvar,  rhsshift,  rhsord,  sistshift,  sisteord; 

u32  enmask; 

struct  eqSymbol  *eq; 

u8  eqendret,  lokrv,  rv=l; 

maxvar=largestSetBit (lc, NWORDS) ; 
enmask=  ( 1«  (maxvar&Oxlf )  )  ; 
enord=maxvar>>5 ; 

for(eq=S->E;  eq<S->E+S->neq;  ++eq) { 
eqendret=0 ; 

for(i=0;  i<eq->nlin;  ++i){ 

if(eq->A[i] [enord] Senmask) { //rekke  i  i  eq  har  variabelen  som  skal 
erstattes 

eqendret=l; 

for ( j  =0 ;  j  <NWORDS ;  ++ j ) 
eq->A[i]  [ j ] A=lc  [  j ] ; 
if (rhv) { //maa  endre  alle  rhs 
rhsshift=l<< (i&Oxlf ) ; 
rhsord=i»5 ; 

for(j=0;  j <eq->originalNRHS ;  ++ j ) { 
if (eq->RHSexists [ j ] ) 

eq->b[j] [rhsord] A=rhsshift; 

} 

} 

} 

} 

if (eqendret) { 

lokrv=trimEquation (eq)  ; 
if (lokrv==0) { 

printf ( "trimeq  returnerer  0  for  ligning  %d\n" , eq->eqnr ) ; 
return  0; 

} 

if ( lokrv==2 ) 
rv=2  ; 

} 

} 

depositLinComb (S->linbank, lc,  rhv)  ; 

//substituted  equation  added  to  the  bank 

return  rv; 


void  trimSystem ( struct  system  *S) { 

/*  Removes  equations  with  no  information  in  them.  */ 
struct  eqSymbol  *eq; 
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for(eq=S->E;  eq<S->E+S->neq;  ++eq) { 

if (eq->nlin<32  &&  ( l<<eq->nlin) ==eq->nrhs ) { //ligning  uten 
inf ormas j  on 

if (eq<S->E+S->neq-l ) 

*eq=S->E [S->neq-l]  ; 
eq— ; 

S->neq-- ; 

} 

} 

} 

void  initLinBank ( struct  eqSymbol  *leq) { 

leq->A=(u32  ** ) malloc (2 *NVAR* sizeof (u32  *)); 

leq->b=(u32  **) malloc ( sizeof (u32  *)); 

leq->b [0] = (u32  *) calloc (2*NWORDS, sizeof (u32) ) ; 

leq->nlin=0 ; 

leq->nrhs=l ; 

leq->nw=0 ; 

leq->nn=0 ; 

leq->eqnr=-l  ; 

leq->originalNRHS=l  ; 

leq->RHSexists= (u8  * ) malloc (sizeof (u8 ) )  ; 
leq->RHSexists [  0 ] =1 ; 


u8  solved ( struct  system  *S) { 
int  i; 


if ( S->neq==0 ) 
return  1; 

if(S->neq==l  &&  S->E  [ 0 ]  . nrhs>=l ) 
return  1; 

for(i=0;  i<S->neq;  ++i){ 
if  (S->E [ i ]  ,nrhs!=l) 
return  0; 

} 

return  1; 


void  solveLinSystem ( struct  eqSymbol  *e)  { 
int  i,  j,  k,  sb,  minb,  maxb,  r; 
u32  *tmp; 
u8  ibit,  rbit; 

for(i=0;  i<e->nlin;  ++i){ 
minb=NVAR; 

for(j=i;  j<e->nlin;  ++ j ) { 

sb=smallestSetBit (e->A[ j ] ,NW0RDS) ; 
if(sb>=0  &&  sbkminb) { 
minb=sb; 
r=j  ; 

} 

if (minb==i ) 
j=e->nlin; 

} 
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if  (e->b  [  0  ]  [i»5]  &  (1«  (i&Oxlf)  )  ) 
ibit=l ; 
else 

ibit=0 ; 

if ( r ! =i ) { / /trenger  a  bytte  om  pa  rader 
tmp=e->A [i ] ; 
e - >A [ i ] =e - >A [ r ] ; 
e->A [r] =tmp; 

if  (e->b  [ 0 ]  [r»5]  &  (1«  (r&Oxlf)  )  ) 
rbit=l ; 
else 
rbit=0 ; 

if (ibitArbit) { 

e->b  [  0  ]  [i»5]  A=  (1«  (i&Oxlf)  )  ; 
ibit A=1 ; 

e->b  [  0  ]  [r»5]  A=  (1«  (r&Oxlf)  )  ; 

} 

} 

for(j=i+l;  j<e->nlin;  ++j) {//lager  0  under  ledende  1 ' er 
if (e->A[ j ] [minb>>5] & (1<< (minb&Oxlf ) ) ) { 
f or ( k=minb>>5 ;  k<NWORDS;  ++k) 
e - > A [ j ] [k] A=e->A[i] [k] ; 
if (ibit) 

e->b[0]  [j>>5]A=(l<<(j&(0xlf) ) ) ; 

} 

} 

} //triangularisert,  starter  tilbakesubstitus j on 

for (i=e->nlin-l;  i>0;  — i ) { 

maxb=largestSetBit (e->A [ i ] , NWORDS )  ; 
if  (e->b  [  0  ]  [i»5]  &  (1«  (i&Oxlf)  )  ) 
ibit=l ; 
else 

ibit=0 ; 

for ( j  =i  —  1 ;  j  >=0 ;  --j){ 

if  (e->A[  j  ]  [maxb>>5]  &  (1«  (maxb&Oxlf )  )  )  { 
for ( k=0 ;  k<NWORDS;  ++k) 
e - > A [ j ] [k] A=e->A[i] [k] ; 
if (ibit) 

e->b  [  0  ]  [j»5]  A=(l«(j&0xlf)  )  ; 

} 

} 

} 

} 

void  copySystem ( struct  system  *S,  struct  system  *kopiS) { 
int  i,  j,  k; 

kopiS->neq=S->neq; 

kopiS->E= ( struct  eqSymbol  *) malloc ( kopiS->neq* sizeof ( struct 
eqSymbol) ) ; 

kopiS->linbank= ( struct  eqSymbol  *) malloc (sizeof (struct  eqSymbol)) 
initLinBank ( kopiS->linbank) ; 
for(i=0;  i<S->linbank->nlin;  ++i){ 
for ( j  =0 ;  j  <NWORDS ;  ++ j  ) 
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kopiS->linbank->A [ i ] [ j ] =S->linbank->A [ i ] [ j ] ; 

} 

for(i=0;  i<S->linbank->nw;  ++i) 

kopiS->linbank->b [ 0 ]  [ i ] =S->linbank->b [ 0 ]  [i]  ; 
kopi S - >1 inban k->nlin=S->l inban k->nl in  ; 
kopi S->1 inban k->nw=S->l inban k->nw; 
for(i=0;  i<S->neq;  ++i){ 

kopiS->E[i] . nlin=S->E [ i ] .nlin; 
kopiS->E[i] . nrhs=S->E [ i ] .nrhs; 

kopiS->E[i] . originalNRHS=S->E [i] .originalNRHS; 
kopiS->E[i] ,nw=S->E[i] .nw; 
kopiS->E[i] .nn=0; 

kopiS->E [i] . A= (u32  **) malloc (kopiS->E [i]  ,nlin*sizeof  (u32  *)); 
for(j=0;  j <kopiS->E [ i ]  . nlin;  ++ j  )  { 

kopiS->E [i].A[j]  =  (u32  *) malloc (NWORDS*sizeof (u32) )  ; 

for ( k=0 ;  k<NWORDS;  ++k) 

kopiS->E [i]  . A [ j ]  [ k] =S— >E [ i ]  .A[ j ]  [k]  ; 

} 

kopiS->E [ i ] . b= (u32  **) malloc ( kopiS->E [ i ]. nrhs*sizeof (u32  *)); 
for(j=0;  j <kopiS->E [ i ]  . nrhs ;  ++ j  )  { 

kopiS->E [i].b[j]=(u32  * ) malloc ( kopiS->E [ i ] . nw*sizeof (u32 ) ) ; 
for(k=0;  k<kopiS->E [ i ] . nw;  ++k) 
kopiS->E [ i ] ,b[j] [ k] =S->E [ i ] -b [ j ] [k] ; 

} 

kopiS->E[i] . RHSexists= (u8  *) malloc ( kopiS- 
>E [i] . originalNRHS* sizeof (u8) ) ; 

for(j=0;  j <kopiS->E [ i ]. originalNRHS ;  ++ j ) 

kopiS->E [i] .RHSexists[j]=S->E[i] . RHSexists [ j ] ; 
kopiS->E[i] . delSinceExtract=S->E [ i ] . delSinceExtract; 

} 

} 

double  log2sum (double  *T,  int  n)  { 

/*  computes  log2 ( \sum_{ i=0 } A { n-1 }  (2 A { T  [ i ] } ) )  .  Works  also  when  the 
sum  is  greater  than  2A{32}.  */ 
double  suma,  sumb; 
int  na; 

if (n==l) 

return  T  [0] ; 
na=n/ 2 ; 

suma=log2  sum ( T , na ) ; 
sumb=log2sum (T+na, n-na) ; 
if (sumaksumb) { 

if ( ( sumb-2 8 . 0 ) >suma) //only  sumb  contributes 
return  sumb; 
else//compute  exactly 

return  suma+log2 (pow (2, sumb-suma) +1.0) ; 

} 

else{//suma  largest 

if ( ( suma-2 8 . 0 ) >sumb) //only  suma  contributes 
return  suma; 
else//compute  exactly 

return  sumb+log2 (pow (2, suma-sumb) +1.0) ; 
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} 


u8  evaluate (u32  *lc,  u8  *v) { 

/*  evaluates  the  linear  combination  lc  with  the  variables  in  v  and 
returns  its  sum  */ 
int  i; 
u8  sum=0; 


for(i=0;  i<NVAR;  ++i){ 

if  (lc  [i»5]  &  (1«  (i&Oxlf )  )  &&  v[i]) 

sumA=l ; 

} 

return  sum; 

} 


2.  AGREEING. h  (by  Havard  Raddum  and  Igor  Semaev) 


u8  checkLink ( struct  linkSymbol  * 1 S ) { 

/*Checks  if  the  two  equations  linked  in  this  link  agree.  Returns  0 
if  equations  are  inconsistent, 

1  if  equations  agree  and  2  if  equations  disagreed  */ 
u8  rv  - 1 ,  lokrv; 

int  i,  j,  eqdel,  linklndex,  eqOni,  eqlni; 
struct  eqSymbol  *eq,  *eqO,  *eql; 

if ( (log2 (lS->ncells) < (log2 (lS->nlist [0] ->nrhs) /2 ) )  && 

( log2 ( lS->ncells ) < ( log2 (lS->nlist [1 ] ->nrhs ) /2) ) ) {//do  not  expect 
deletions 

for(i=0;  i<lS->ncells ;  ++i) { 
if ( lS->cellExists [i]  )  { 
eqdel=2 ; 

if ( lS->nCover [ 0 ]  [ i ] == 0  &&  lS->nCover [ 1 ]  [i]>0) 
eqdel=l ; 

if ( lS->nCover [ 0 ]  [ i ] >0  &&  lS->nCover [ 1 ]  [ i ] == 0 ) 
eqdel=0 ; 
if (eqdel<2 ) { 
rv=2  ; 

eq=lS->nlist [eqdel]  ; 
for(j=0;  j<eq->nn;  ++ j ) { 
if  (eq->link  [  j  ]  =1S)  { 
linklndex= j ; 
j  =eq->nn+l ; 

} 

} 

if ( j  ==eq->nn)  { 

printf("Fant  ikke  link!?  eqdel=%d,  i=%d\n" , eqdel, i) ; 
printf ( "Linkadresse :  %x,  eq . naboadresse :  %x\n",lS,eq- 

>link [ 0 ]  )  ; 

} 

j=0; 
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while ( lS->nCover [eqdel ] [i]>0) { 
while ( j <eq->originalNRHS  &&  eq- 
>coverID [ linklndex]  [ j ]!  =  (unsigned) i )  //  !!  added  cast 

++j  ; 

if ( j ==eq->originalNRHS ) printf ( " f ant  ikke  RHS  som  dekker  celle 
%x,  13->nCover=%d\n" , i , lS->nCover [eqdel ] [i]); 
lokrv=deleteRHS (eq,  j )  ; 
if ( lokrv==0 ) 
return  0; 

++j  ; 

} 

lS->cellExists [ i ] =0 ; 

} 

} 

} 

} 

else { //expects  deletions,  goes  through  RHS ' s  of  linked  equations 
eqO=lS->nlist [0]  ; 
eql=lS->nlist[l]  ; 
for(i=0;  i<eq0->nn;  ++i){ 
if (eq0->link [ i ] ==1S )  { 
eq0ni=i ; 
i=eq0->nn+l  ; 

} 

} 

if (i==eq0->nn) printf (" (checkLink) Fant  ikke  link  fra  nabo  0\n"); 
for(i=0;  i<eql->nn;  ++i){ 
if (eql->link [i] ==1S)  { 
eqlni=i ; 
i=eql->nn+l ; 

} 

} 

if (i==eql->nn) printf (" (checkLink) Fant  ikke  link  fra  nabo  l\n"); 
for(i=0;  i<eql->originalNRHS;  ++i){ 
if (eql->RHSexists [i]  )  { 

if ( lS->nCover [ 1 ] [eql->coverID [eqlni] [ i ] ] ==0 ) 

printf ("RHS  %d  i  eql  dekker  en  celle  som  link  sier  ikke  blir 
dekket  av  eql\n",i); 

if (lS->nCover [0] [eql->coverID [eqlni] [i] ] ==0) { 
lokrv=deleteRHS (eql,  i)  ; 
if ( lokrv==0 ) 
return  0; 
rv=2 ; 

} 

} 

} 

for(i=0;  i<eqO->originalNRHS ;  ++i){ 
if (eq0->RHSexists [i] ) { 

if ( lS->nCover [ 0 ] [eq0->coverID [eqOni ] [ i ] ] ==0 ) 

printf ("RHS  %d  i  eqO  dekker  en  celle  som  link  sier  ikke  blir 
dekket  av  eq0\n",i); 

if ( lS->nCover [ 1 ] [eq0->coverID [eqOni ] [ i ] ] ==0 ) { 
lokrv=deleteRHS (eqO, i) ; 
if ( lokrv==0 ) 
return  0; 
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rv=2  ; 

} 

} 

} 

} 

return  rv; 

} 

u8  agreeSystem ( struct  system  *S) { 

/*  Agrees  the  whole  system.  Returns  0  if  system  is  inconsistent, 
1  if  system  already  agreed,  or  2  if  deletions  of  RHS ' s  have 
occurred  */ 
int  i; 

u8  rv=l,  lokrv,  changed=l; 
struct  linkSymbol  *lsym; 

while (changed) { 
changed=0 ; 

f or ( lsym=S->L;  lsym<S->L+S->nlink;  ++lsym) { 
lokrv=checkLink ( lsym) ; 
if ( lokrv==2 ) 
changed=l ; 
if (lokrv— =0 )  { 

printf ( "inconsistency  in  link  %d\n",i); 
return  0; 

} 

} 

} 

return  rv; 

} 


} 


3.  GLUING. h  (by  Havard  Raddum  and  Igor  Semaev) 

/*  Methods  inmplementing  gluing  of  equations,  before  linking  */ 

u8  maskedRHSlEqualRHS2 (u32  *rhsl,  u32  *rhs2,  int  nw,  u32  *maske) { 
int  t; 

t=nw-l ; 

while (t>=0  &&  ( (rhsl [t] Smaske [t] ) == (rhs2 [t] Smaske [t] ) ) ) 
t--; 
if (t<0) 
return  1; 
else 

return  0; 


u8  maskedRHSlBiggerThanRHS2 (u32  *rhsl,  u32  *rhs2,  int  nw,  u32  *maske) { 
int  t; 

t=nw-l ; 

while (t>=0  &&  (rhsl [t] Smaske [t] )== (rhs2 [t] Smaske [t] ) ) 
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t--; 

if (t>=0  &&  (rhsl [t] Smaske [t] ) > (rhs2 [t] Smaske [t] ) ) 
return  1; 
else 

return  0; 


void  mergeSortMaskedRightHandSides (u32  **RHS,  int  nicl,  int  nw,  u32 
*maske) { 

int  niacl,  nibcl,  i,  j,  t; 
u32  *tmp,  **aRHS,  **bRHS,  **dl; 

if (nicl==2  &&  mas kedRHSl Bigger ThanRHS 2 (RHS [ 0 ] , RHS [ 1 ] , nw, maske) ) { 
tmp=RHS [ 0 ] ; 

RHS [0] =RHS [1] ; 

RHS [1] =tmp; 

} 

if (nicl>2 ) { 

niacl=nicl/ 2 ; 
nibcl=ni cl -niacl ; 
aRHS=RHS ; 
bRHS=RHS+niacl; 

mergeSortMaskedRightHandSides (aRHS, niacl, nw, maske) ; 
mergeSortMaskedRightHandSides (bRHS , nibcl , nw, maske) ; 
dl=(u32  ** ) malloc (nicl* sizeof (u32  *)  )  ; 
i=j=t=0; 

while ( ikniacl  &&  j <nibcl ){ //ingen  er  ferdige,  ma  sammenligne 
if (mas kedRHSl Bigger ThanRHS 2 (aRHS [ i ] , bRHS [ j ] , nw, maske) ) 
dl [ t++ ] =bRHS [ j  ++ ] ; 
else 

dl [ t++ ] =aRHS [ i++ ] ; 

} 

if ( i==niacl ) { //aRHS  ble  ferdig  f0rst,  kopierer  resten  av  bRHS 
for(i=j;  iCnibcl;  ++i) 
dl [t++] =bRHS [i] ; 

} 

else {//bRHS  ble  ferdig  fprst,  kopierer  resten  av  aRHS 
for(j=i;  jcniacl;  ++ j ) 
dl [t++] =aRHS [ j ] ; 

} 

//kopierer  tilbake  i  RHS 
for(i=0;  i<nicl;  ++i) 

RHS [i]=dl [i]  ; 
free (dl ) ; 

} 

} 

void  computeMaskedSortedRHS ( struct  eqSymbol  *eql,  struct  eqSymbol  *eq2, 

u32  **Ubl,  u32  **Ub2,  u32  *maske,  u32  **U,  int 

ncommon) { 

/*  Expands  RHS ' s  in  eql  and  eq2  so  they  can  be  glued  together.  maske 
shows  which  bits  that  must  be  equal  when  gluing.  */ 
u32  *lokb; 

int  i,  j,  t,  sumnl,  sumnlw,  skev; 
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sumnl=eql->nlin+eq2->nlin; 
sumnlw= (sumnl  +  31)  >>5; 

lokb= (u32  * ) malloc ( sumnlw*sizeof (u32 )  )  ; 
t=0; 

for(i=0;  i<eql->originalNRHS;  ++i ){ //computes  Ub  for  all  right-hand 
sides  b  in  El 

if (eql->RHSexists  [i] )  { 

for(j=0;  j<eql->nw;  ++ j ) 
lokb [ j ] =eql->b [i]  [  j  ]  ; 
for ( j=eql->nw;  j<sumnlw;  ++ j ) 
lokb [ j ] =0 ; 

Uxb (U,  sumnl,  sumnl,  lokb,  Ubl [t++] ) ; 

} 

} 

if (t ! =eql->nrhs) printf ( "eql->nrhs=%d,  t=%d\n" , eql->nrhs, t) ; 
t=0; 

for(i=0;  i<eq2->originalNRHS;  ++i ){ //computes  Ub  for  all  right-hand 
sides  b  in  E2 

if (eq2->RHSexists  [i] )  { 

for(j=0;  j<eql->nw;  ++ j ) 
lokb [ j ] =0 ; 

skev= (eql->nlin) &0xlf ; 
if ( skev==0 ) { 

for ( j=eql->nw;  j<sumnlw;  ++ j ) 
lokb [ j ] =eq2->b [i] [j-eql->nw]; 

} 

else  { 

for(j=0;  j<eq2->nw;  ++ j ) { 

lokb  [eql->nw+j  -1]  |  =  (eq2->b  [i]  [  j  ]  «skev)  ; 
if (sumnlw> ( j+eql->nw) ) 

lokb [eql->nw+j ] =eq2->b [i] [ j ] >> (32-skev) ; 

} 

} 

Uxb (U, sumnl, sumnl, lokb, Ub2 [t++] ) ; 

} 

} 

if (t ! =eq2->nrhs) printf ( "eq2->nrhs=%d,  t=%d\n" , eq2->nrhs, t) ; 

free (lokb) ; 

i=sumnl-ncommon; 

j=0; 

while ( i>=32 ) { 
maske [ j++] =0; 
i-=32 ; 

} 

if ( j  <sumnlw)  { 

maske [j]=(0xffffffff<<i) ; 
for(i=j+l;  i<sumnlw;  ++i) 
maske [i]=0xffffffff; 
if (sumnl&Oxlf) 

maske [sumnlw-1] &= (1<< (sumnl&Oxlf) ) -1; 

} 

//mask  in  place 
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mergeSortMaskedRightHandSides (Ubl, eql->nrhs, sumnlw, maske) ; 
merge SortMaskedRightHandS ides (Ub2 , eq2->nrhs , sumnlw, maske) ; 


int  nRHSwhenGlued ( struct  eqSymbol  *eql,  struct  eqSymbol  *eq2){ 

/*  Computes  number  of  RHS  if  gluing  eql  and  eq2  */ 
u32  **Ubl,  *  *Ub2 ,  *maske,  **U; 

int  i,  j,  nn=0,  bpl,  ppl,  bp2,  pp2,  prod,  nw,  nl,  nc; 

nl=eql->nlin+eq2->nlin; 
nw= (nl+31) >>5; 

U=(u32  *  * ) malloc (nl* sizeof (u32  *)); 
nc=computeU (eql , eq2 , U) ; 
if (nc==0 ) { 

for(i=0;  i<nl;  ++i) 
free (U [i] )  ; 
free (U) ; 

if ( (MAXINT/eql->nrhs) <eq2->nrhs) 
return  MAXINT; 
else 

return  eql->nrhs*eq2->nrhs; 

} 

Ubl=(u32  **) malloc (eql->nrhs* sizeof (u32  *)); 
for(i=0;  i<eql->nrhs;  ++i) 

Ubl [i] = (u32  *)calloc(nw,sizeof(u32) ) ; 

Ub2=(u32  **) malloc (eq2->nrhs*sizeof (u32  *)); 
for(i=0;  i<eq2->nrhs;  ++i) 

Ub2 [i]=(u32  *)calloc(nw,sizeof(u32) ) ; 
maske= (u32  * ) malloc (nw* sizeof (u32 ) ) ; 

computeMaskedSortedRHS (eql , eq2 , Ubl , Ub2 , maske,  U,  nc) ; 

for(i=0;  i<nl;  ++i) 
free (U  [i] )  ; 
free (U) ; 

bp 1 =pp 1 =bp2 =pp2 = 0 ; 

while (bpl<eql->nrhs  &&  bp2<eq2->nrhs  &&  nn<MAXINT) { 
while (bpl<eql->nrhs  &&  bp2<eq2->nrhs  && 

!maskedRHSlEqualRHS2 (Ubl [bpl] , Ub2 [bp2] , nw, maske) ) { 

if (maskedRHSlBiggerThanRHS2 (Ubl [bpl ] , Ub2 [bp2 ] , nw, maske) ) 

bp2++; 

else 

bpl++; 

}//her  er  maskert  Ubl [bpl]  og  Ub2[bp2]  like,  eller  en  av  listene  er 
nadd  til  endes 
pp 1 =bp 1 ; 
pp2=bp2 ; 

while (ppl<eql->nrhs  && 

maskedRHSlEqualRHS2 (Ubl [bpl ] , Ubl [ppl ]  ,  nw,  maske) ) 
ppl++; 

while (pp2<eq2->nrhs  && 

maskedRHSlEqualRHS2 (Ub2 [bp2 ]  ,  Ub2 [pp2 ]  ,  nw,  maske) ) 

pp2++; 

//Blokkene  fra  bpl  og  bp2  til  ppl  og  pp2  er  identiske  maskert 
if ( (ppl-bpl ) >0  &&  MAXINT/ (ppl-bpl) < (pp2-bp2) ) 
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nn=MAXINT; 

prod= (ppl-bpl ) * (pp2-bp2 ) ; 
if (nn>MAXINT-prod) 
nn=MAXINT; 
if (nn<MAXINT) 
nn+=prod; 
bp 1 =pp 1 ; 
bp2=pp2 ; 

} 

for(i=0;  i<eql->nrhs;  ++i) 
free (Ubl  [  i ] )  ; 
free (Ubl ) ; 

for(i=0;  i<eq2->nrhs;  ++i) 
free (Ub2  [  i ] )  ; 
free (Ub2 ) ; 
free (maske) ; 

return  nn; 

} 

int  estimateNRHSwhenGlued ( struct  eqSymbol  *eql,  struct  eqSymbol  *eq2){ 
int  i,  enn,  nc,  suranl; 
u32  **M; 

sumnl=eql->nlin+eq2->nlin ; 

M=(u32  ** ) malloc ( sumnl*sizeof (u32  *)); 
for(i=0;  i<eql->nlin;  ++i) 

M  [  i ] =eql->A [ i ] ; 

for (i=eql->nlin;  i<sumnl;  ++i) 

M[i] =eq2->A[i-eql->nlin]  ; 
nc=sumnl-rank (M, sumnl , NVAR) ; 
free (M) ; 

if ( (log2 (eql->nrhs) +log2 (eq2->nrhs) - (double) nc) <30 . 9) { 
if ( (log2 (eql->nrhs) +log2 (eq2->nrhs) - (double) nc) >0 . 0)  { 
if (eq2->nrhs>eql->nrhs) 

enn= (eql->nrhs>> (nc/ 2 ) ) * (eq2->nrhs>> (nc- (nc/ 2 ) ) ) ; 
else 

enn= (eq2->nrhs>> (nc/ 2 ) ) * (eql->nrhs>> (nc- (nc/ 2 ) ) ) ; 
if (enn==0) 
enn=l ; 

} 

else 
enn=l ; 

} 

else 

enn=MAXINT; 
return  enn; 

} 

void  glue (struct  eqSymbol  *eql,  struct  eqSymbol  *eq2,  struct  eqSymbol 
*geq) { 

/*  Glues  together  eql  and  eq2  into  geq,  eql  and  eq2  should  not  be 
used  afterwards.  */ 

int  i,  j,  k,  bpl,  ppl,  bp2,  pp2,  ncommon,  t,  estnnew; 
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u32  *maske,  **Ubl,  **Ub2,  **U,  **nyA; 

Ubl=(u32  *  * ) malloc (eql->nrhs* sizeof (u32  *)); 

Ub2=(u32  **) malloc (eq2->nrhs*sizeof (u32  *)); 

geq->nlin=eql->nlin+eq2->nlin; 

geq->nw= (geq->nlin+31 ) >>5; 

geq->A=(u32  **) malloc (geq->nlin*sizeof (u32  *)); 
for (1=0;  i<eql->nlin;  ++i) 
geq->A[i] =eql->A[i] ; 
for(i=0;  i<eq2->nlin;  ++i) 

geq->A[eql->nlin+i] =eq2->A[i]  ; 

nyA=(u32  **) malloc (geq->nlin* sizeof (u32  *)); 

U=(u32  **)  malloc (geq->nlin*sizeof (u32  *)); 
for(i=0;  i<geq->nlin;  ++i) 

nyA[i]=(u32  *) calloc (NWORDS, sizeof (u32) ) ; 
for(i=0;  i<eql->nrhs;  ++i) 

Ubl [i] = (u32  *)calloc (geq->nw, sizeof (u32) ) ; 
for(i=0;  i<eq2->nrhs;  ++i) 

Ub2 [i]=(u32  *) calloc (geq->nw, sizeof (u32) ) ; 
maske= (u32  * ) malloc (geq->nw* sizeof (u32 ) ) ; 

ncommon=computeU (eql , eq2  ,  U) ; 

UAX (U, geq->A, nyA, geq->nlin, NWORDS ) ; 

free (geq->A) ; 

geq->A=nyA; 

computeMaskedSortedRHS (eql , eq2 , Ubl , Ub2 , maske, U, ncommon) ; 
estnnew= (eql->nrhs>> (ncommon/ 2 ) ) * (eq2->nrhs>> (ncommon-ncommon/ 2 ) ) ; 
if (estnnewkl ) 
estnnew=l ; 

//printf ( "estnnew=2  A ( % . 3f ) ,  ", log2 (estnnew) ) ; 

geq->b=(u32  **) malloc (estnnew*sizeof (u32  *)); 
bpl=ppl=bp2=pp2=t=0  ; 

while (bpl<eql->nrhs  &&  bp2<eq2->nrhs) { 
while (bpl<eql->nrhs  &&  bp2<eq2->nrhs  && 

!maskedRHSlEqualRHS2 (Ubl [bpl] , Ub2 [bp2] , geq->nw, maske) ) { 

if (maskedRHSlBiggerThanRHS2 (Ubl [bpl ] , Ub2 [bp2 ] , geq->nw, maske) ) 

bp2++; 

else 

bpl++; 

}//her  er  maskert  Ubl [bpl]  og  Ub2[bp2]  like,  eller  en  av  listene  er 
nadd  til  endes 
pp 1 =bp 1 ; 
pp2=bp2 ; 

while (ppl<eql->nrhs  &&  maskedRHSlEqualRHS2 (Ubl [bpl ] , Ubl [ppl ] , geq- 
>nw, maske) ) 
ppl++; 

while (pp2<eq2->nrhs  &&  maskedRHSlEqualRHS2 (Ub2 [bp2 ] , Ub2 [pp2 ] , geq- 
>nw, maske) ) 

pp2++; 

//Blokkene  fra  bpl  og  bp2  til  ppl  og  pp2  er  identiske  maskert, 
limer  alle  par 

for(i=bpl;  i<ppl;  ++i) { 
for ( j  =bp2 ;  j<pp2;  ++ j )  { 
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geq->b [t] = (u32  * ) malloc (geq->nw*sizeof (u32 ) ) ; 
for(k=0;  k<geq->nw;  ++k) 

geq->b [t] [k] =Ubl [i] [k] AUb2 [ j ] [k] ; / /faktisk  liming 
t++; 

if (t==estnnew) { //needs  more  memory  and  a  better  estimate 
if (i==0) 

estnnew*=2 ; 
else  { 

if (i==eql->nrhs) 
estnnew+= (pp2- j ) ; 
else 

estnnew= (int) ( (double) estnnew/ ( (double) i/ (double) eql- 

>nrhs) ) ; 

} 

geq->b=(u32  **) realloc (geq->b, estnnew* sizeof (u32  *)); 

} 

} 

} 

bpl=ppl; 
bp2=pp2 ; 

} 

geq->originalNRHS=t  ; 
geq->nrhs=t; 

geq->b= (u32  **) realloc (geq->b, geq->originalNRHS* sizeof (u32  *)); 
geq->RHSexists= (u8  * ) malloc (geq->originalNRHS* sizeof (u8 ) ) ; 
for(i=0;  i<geq->originalNRHS ;  ++i) 
geq->RHSexists [ i J  -1 ; 

if (ncommon>0 ){ //trimmer  bort  lineasrt  avhengig  informasjon 
for (i=geq->nlin-ncommon;  i<geq->nlin;  ++i) 
free (geq->A [ i ] ) ; 
geq->nlin-=ncommon; 

geq->A=(u32  **) realloc (geq->A, geq->nlin*sizeof (u32  *)); 
geq->nw= (geq->nlin+31 )  >>5; 
if ( ( (geq->nlin+31 ) >>5) <geq->nw) { 
geq->nw= (geq->nlin+31 ) >>5; 
for(i=0;  i<geq->nrhs;  ++i) 

geq->b [ i ] = (u32  * ) realloc (geq->b [ i ] , geq->nw* sizeof (u32 ) ) ; 

} 

} 

geq->nn=0 ; 
free (maske) ; 

for(i=0;  i<eql->nrhs;  ++i) 
free (Ubl [ i ] ) ; 
free (Ubl ) ; 

for(i=0;  i<eq2->nrhs;  ++i) 
free (Ub2 [ i ] ) ; 
free (Ub2 ) ; 


u8  packSystem ( struct  system  *S,  int  th) { 

/*  Tries  to  glue  together  as  many  equations  as  possible,  with  the 
restriction  that  the  number  of 

right-hand  sides  in  the  glued  equation  should  not  be  estimated  to 
be  larger  than  th. 
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Returns  0  if  inconsistencies  are  found,  1  otherwise.  */ 
struct  eqSymbol  *nyE,  *tmp; 

int  i,  j,  nnew,  minnew,  mini,  min j ,  nyneq=0,  **gluetab; 
u8  *limt,  tryglue,  trypack; 

nyE= (struct  eqSymbol  *) malloc ( S->neq* sizeof ( struct  eqSymbol)); 
tmp= (struct  eqSymbol  *) malloc ( sizeof ( struct  eqSymbol)); 
limt= (u8  *)calloc (S->neq, sizeof (u8) ) ; 

gluetab=(int  **) malloc ( S->neq*sizeof ( int  *)); 
for(i=0;  i<S->neq;  ++i) 

gluetab[i]=(int  *)calloc (S->neq, sizeof (int) ) ; 
for(i=0;  i<S->neq;  ++i){ 
for(j=i  +  l;  j<S->neq;  ++ j  )  { 

glue tab [ i ] [ j ] =estimateNRHSwhenGlued ( S->E+i , S->E+ j ) ; 
if (glue tab [ i ] [ j ] ==0 ) { 
j  =S->neq; 
i=S->neq+l ; 

} 

} 

} 

if ( i==S->neq) 
tryglue=l ; 

else { / /funnet  inkonsistens 
tryglue=2 ; 

printf  ("inconsistency  foundin'')  ; 

} 

while (tryglue==l ) { 
minnew=th+l ; 

for(i=0;  i<S->neq;  ++i){ 
if ( ! limt [ i ] ) { 

for(j=i  +  l;  j<S->neq;  ++ j  )  { 
if ( ! limt [ j ] ) { 

nnew=gluetab [ i ] [j]; 
if (nnew<minnew) { 
minnew=nnew; 
mini=i ; 
min j  = j ; 

} 

} 

} 

} 

}//her  er  billigste  liming  funnet,  hvis  mulig  aa  lime  under  th 
if (minnew<=th) { //mulig  aa  lime 

glue ( S->E+mini , S->E+min j  ,  nyE+nyneq)  ; 

nyE[nyneq] ,eqnr=nyneq; 

limt [mini ] =1  ; 

limt [min j ] =1  ; 

trypack=l ; 

while (trypack) { 

minnew=th+l ; 

for(i=0;  i<S->neq;  ++i){ 
if ( ! limt [ i ] ) { 

nnew=estimateNRHSwhenGlued (S->E+i,  nyE+nyneq)  ; 
if (nnew==0 ) { 
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tryglue=2 ; 

printf ("Inconsistency  Found\n") ; 
i=S->neq; 

} 

else  { 

if (nnew<minnew) { 
minnew=nnew; 
mini=i ; 

} 

} 

} 

} 

if (tryglue==l  &&  minnew<=th) { //mulig  aa  putte  paa  en  til 
//printf ("  og  %d\n",mini); 
glue (nyE+nyneq, S->E+mini, tmp) ; 
limt [mini ] =1 ; 

deleteEquation (nyE+nyneq) ; 
nyE [nyneq] =*tmp; 

} 

else 

trypack=0 ; 

} 

nyneq++; 

} 

else 

tryglue=0 ; 

}//har  pakket  saa  godt  som  mulig,  kopierer  ulimte  symboler 
if (tryglue<2 ) { 

for(i=0;  i<S->neq;  ++i){ 
if (limt [i] ) 

deleteEquation ( S->E+i ) ; 
else  { 

nyE [nyneq] =S->E [ i ]  ; 
nyE [nyneq] . eqnr=nyneq; 
nyneq++; 

} 

} 

free (S->E) ; 

S->E=nyE; 

S->neq=nyneq; 
for(i=0;  i<S->neq;  ++i) 

S->E[i] ,eqnr=i; 

} 

free (limt) ; 

for(i=0;  i<S->neq;  ++i) 
free (gluetab [ i ] ) ; 
free (gluetab) ; 
if (tryglue==2) 
return  0; 
return  1; 
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4.  LINKING. h 


/*  Methods  for  linking  equations  in  system  */ 

void  zeroPrepend (u32  *b,  int  lb,  int  nzeros,  u32  *x) { 
int  xw,  nOw,  lbw,  skev,  i; 

nOw=  (nzeros»5)  ; 
xw= ( lb+nzeros  +  31 ) >>5 ; 
lbw= (lb+31) >>5; 
for(i=0;  i<xw;  ++i) 
x  [  i  ]  =  0  ; 

skev= (nzeros&Oxlf )  ; 
if (skev) { 

for(i=0;  i<lbw;  ++i){ 

x  [  i+nOw]  |  =  (b  [  i  ]  «skev)  ; 

if ( (i+nOw) < (xw-1) ) 

x [i+nOw+1] | = (b [i] >> (32 -skev) ) ; 

} 

} 

else  { 

for(i=0;  i<lbw;  ++i) 
x [ i+nOw] =b [ i ] ; 

} 

} 

void  zeroAppend (u32  *b,  int  lb,  int  nzeros,  u32  *x) { 
int  i,  xw,  lbw; 

lbw= (lb+31) >>5; 
xw= (lb+nzeros+31) >>5; 
for(i=0;  i<lbw;  ++i) 
x [ i ] =b  [  i ] ; 

for(i=lbw;  i<xw;  ++i) 

x [ i ] =0 ; 

} 

int  collapseCelllDlist (struct  bitVector  *vl,  int  nil){ 
int  bp=0,  pp=l; 

while (ppcnil) { 

while (pp<nil  &&  vOEqualVl (vl [bp] , vl [pp] ) ) 
free (vl [pp++] . v) ; 
if (pp<nil ) 

vl [++bp] — v 1 [pp++] ; 

} 

return  bp+1; 


int  findCelllndex (struct  bitVector  *vl,  int  nil,  struct  bitVector 
target)  { 

int  bunn,  topp,  midt; 


if (vOGreaterThanVl (vl [0] , target) ) //target  smaller  than  smallest  in  vl 
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return  -1; 

if (vOGreaterThanVl (target ,  vl [nil-1 ])) //target  greater  than  largest  in 
return  -1; 

if (vOEqualVl (vl [0] , target) ) 
return  0; 

if (vOEqualVl (vl [nil-1 ] , target) ) 
return  nil-1; 

//Know  target  is  not  inside  boundary  of  vl 
bunn=0 ; 
topp=nil-l ; 
while ( (topp-1) >bunn) { 
midt= (bunn+topp) / 2 ; 
if (vOEqualVl (vl [midt] , target) ) 

return  midt; 

if (vOGreaterThanVl (vl [midt] , target) ) 
topp=midt; 
else 

bunn=midt; 

} 

return  -1; 


void  establishSmallLink ( struct  eqSymbol  *eq0,  struct  eqSymbol  *eql, 
struct  linkSymbol  *ls,  u32  **M0,  int  commondim) { 

/*  Spaces  spanned  by  A-matrices  in  equations  overlap  in  few 
dimensions,  constructs  a  link  with  2Ar  cells  */ 
int  i,  j,  sumnl,  sumnlw; 
u32  *x,  *ci; 

sumnl=eqO->nlin+eql->nlin; 
sumnlw= (sumnl+31) >>5; 
ls->nlist[0]=eq0; 
ls->nlist[l] =eql ; 
eq0->link [eq0->nn] =ls; 
eql->link [eql->nn] =ls; 

ls->ncells= ( l<<commondim) ; //commondim  should  be  small  enough  for  this 
ls->cellExists= (u8  * ) malloc ( ls->ncells*sizeof (u8 ) )  ; 
for(i=0;  i<2;  ++i) 

ls->nCover [i] = (int  *)calloc (ls->ncells, sizeof (int) ) ; 
eqO->coverID [eq0->nn]  =  (u32  * ) malloc (eqO->originalNRHS* sizeof (u32 ) )  ; 
eql->coverID [eql->nn]  =  (u32  * ) malloc (eql->originalNRHS* sizeof (u32 ) ) ; 
x— (u32  * ) malloc ( sumnlw* sizeof (u32  )  )  ; 
ci=(u32  *) malloc ( sizeof (u32 ))  ; 

/*  Establish  connections  for  eqO  */ 
for(i=0;  i<eqO->originalNRHS ;  ++i){ 
if (eqO->RHSexists [i]  )  { 

zeroAppend (eq0->b [i] , eq0->nlin, eql->nlin,  x)  ; 

Uxb (MO , commondim, sumnl ,  x,  ci )  ; 
ls->nCover [0] [ci [0] ] ++; 
eqO->coverID [eq0->nn]  [ i ] =ci  [ 0 ] ; 

} 
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} 

/*  Establish  connections  for  eql  */ 
for(i=0;  i<eql->originalNRHS;  ++i){ 
if (eql->RHSexists [i] ) { 

zeroPrepend (eql->b [i] , eql->nlin, eqO->nlin,  x) ; 
Uxb (MO , comraondim, sumnl , x,  ci ) ; 
ls->nCover [1] [ci [0] ] ++; 
eql->coverID [eql->nn] [i] =ci [0] ; 

} 

} 

for(i=0;  i<ls->ncells ;  ++i) { 

if ( ls->nCover [ 0 ]  [ i ] >  0  |  |  ls->nCover [ 1 ]  [i]>0) 
ls->cellExists [ i ] =1 ; 
else 

ls->cellExists [ i ] =0 ; 

} 

eq0->nn++; 

eql->nn++; 

free (x) ; 
free  (ci) ; 


void  establishBigLink ( struct  eqSymbol  *eqmin,  struct  eqSymbol  *eqmax, 
struct  linkSymbol  *ls,  u32  **M0,  int  commondim) { 

/*  Constructs  a  link  where  the  equations  overlap  in  too  many 
dimensions.  The  (number  of)  cells  in  the  link  is  based  on 
the  right-hand  sides  from  the  smallest  equation  */ 
int  i,  j,  t,  sumnl,  sumnlw,  lew,  newNcell,  cellindex; 
u32  *x; 

struct  bitVector  ci; 

sumnl=eqmin->nlin+eqmax->nlin; 
sumnlw= (sumnl+31) >>5; 
lcw= (commondim+31 ) >>5; 
ls->nlist [0] =eqmin; 
ls->nlist [1] =eqmax; 
eqmin->link [eqmin->nn] =ls  ; 
eqmax->link [eqmax->nn] =ls; 
ls->ncells=eqmin->nrhs  ; 

if (ls->ncells==0) {printf (" (establishBigLink) Equations  %d  and  %d 
inconsistent\n" , eqmin->eqnr , eqmax->eqnr ) ; exit ( 0 ) ; } 
ls->cellExists= (u8  * ) malloc ( ls->ncells* sizeof (u8 ) )  ; 
ls->cellID= (struct  bitVector  *) malloc (ls->ncells*sizeof (struct 
bitVector) ) ; 

for(i=0;  i<2;  ++i) 

ls->nCover [i] = (int  *) calloc (ls->ncells, sizeof (int) ) ; 
eqmin->coverID [eqmin->nn] = (u32  * ) malloc (eqmin- 
>originalNRHS*sizeof (u32) ) ; 

eqmax->coverID [eqmax->nn] = (u32  * ) malloc (eqmax- 
>originalNRHS*sizeof (u32) ) ; 

x— (u32  * ) malloc ( sumnlw* sizeof (u32 ) ) ; 
ci . v— (u32  * ) malloc ( lew* sizeof (u32 ) ) ; 
ci . length=commondim; 
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ci . wl=lcw; 

/*  Establish  connections  for  eqmin  */ 
t=0; 

for(i=0;  i<eqmin->originalNRHS ;  ++i){ 
if (eqmin->RHSexists [i]  )  { 

zeroAppend (eqmin->b [i] , eqmin->nlin, eqmax->nlin,  x) ; 

Uxb (MO , commondim, sumnl , x, ci . v) ; 
ls->cellID [t] ,wl=lcw; 
ls->cellID [t] . length=commondim; 

ls->cellID [ t ]  . v— (u32  * ) malloc ( lcw*sizeof (u32 ) )  ; 
for(j=0;  j<lcw;  ++ j ) 
ls->cellID [t] . v [ j ] =ci . v [ j ] ; 
t++; 

} 

}//all  different  celllDs,  with  repetition,  created 

if ( t ! =eqmin->nrhs ) {printf ( " (establishBigLink) More  existing  RHS ' s  (%d) 
in  eqmin  than  eqmin->nrhs=%d  says!\n", 

t , eqmin- >nrhs ) ; exit (0 ) ; } 
mergeSortBitVectors (ls->cellID, ls->ncells)  ; 
newNcell=collapseCellIDlist (ls->cellID, ls->ncells) ; 
if (newNcell<ls->ncells ) { 

ls->ncells=newNcell ; //the  actual  number  of  cells  in  this  link 
ls->cellID= (struct  bitVector  *) realloc (ls->cellID, ls- 
>ncells*sizeof (struct  bitVector) ) ; 

} 

for(i=0;  i<2;  ++i) 

ls->nCover [i] = (int  *)calloc (ls->ncells, sizeof (int) ) ; 
for(i=0;  i<eqmin->originalNRHS ;  ++i){ 
if (eqmin->RHSexists [i]  )  { 

zeroAppend (eqmin->b [i] , eqmin->nlin, eqmax->nlin, x) ; 

Uxb (MO , commondim, sumnl ,  x,  ci  .  v)  ; 

cellindex=f indCelllndex (ls->cellID,  ls->ncells,  ci)  ; 
if (cellindex==-l ) {printf ( " (establishBigLink) fant  ikke 
cell ID ! " ) ; exit  ( 0 ) ;  } 

eqmin- >cover ID [ eqmin- >nn] [ i ] =ce 11 index; 
ls->nCover [ 0 ] [cellindex] ++; 

} 

}//link  to  eqmin  created 
/*  Establish  connections  for  eqmax  */ 
for(i=0;  i<eqmax->originalNRHS ;  ++i){ 
if (eqmax->RHSexists [i] ) { 

zeroPrepend (eqmax->b [i] , eqmax->nlin, eqmin->nlin, x) ; 

Uxb (MO , commondim, sumnl , x, ci . v) ; 

cellindex=f indCelllndex (ls->cellID, ls->ncells,  ci) ; 
if (cellindex==-l ) //RHS  i  does  not  agree  with  eqmin 
deleteRHS (eqmax, i) ; 
else  { 

eqmax- >cover ID [ eqmax- >nn] [ i ] =ce 11 index; 
ls->nCover [ 1 ] [cellindex] ++; 

} 

} 

} //links  created,  no  more  use  for  celllD,  deletes  it 
for(i=0;  i<ls->ncells ;  ++i) 
free (ls->cellID [i] . v) ; 
free (ls->cellID) ; 
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ls->cellExists= (u8  * ) malloc ( ls->ncells*sizeof (u8 ) ) ; 
for(i=0;  i<ls->ncells ;  ++i)  { 

if ( ls->nCover [ 0 ]  [ i ] > 0  |  |  ls->nCover [ 1 ]  [ i ] >0 ) 

ls->cellExists [ i ] =1 ; 
else 

ls->cellExists [ i ] =0 ; 

} 

eqmin->nn++; 

eqmax->nn++; 

free (x) ; 
free (ci . v) ; 


void  tryEstablishLink ( struct  system  *S,  int  eqnrO,  int  eqnrl){ 
struct  eqSymbol  *eqmin,  *eqmax; 
struct  linkSymbol  *ls; 

int  i,  commondim,  suranl; //  ! !  deleted  var:  ,  sumnlw; 
u32  **U; 

if (S->E [eqnrO] . nrhs<S->E [eqnrl ] .nrhs) { 
eqmin=S->E+eqnrO  ; 
eqmax=S->E+eqnrl ; 

} 

else  { 

eqmin=S->E+eqnrl  ; 
eqmax=S->E+eqnrO  ; 

}//eqmin  goes  on  top  when  computing  U! 
sumnl=eqmin->nlin+eqmax->nlin; 

//  !!  deleted  var:  sumnlw= (sumnl+31) >>5; 

U=(u32  **) malloc ( sumnl*sizeof (u32  *)); 
commondim=computeU (eqmin, eqmax,  U) ; 

if (commondim>0) { /*  eqmin  and  eqmax  can  exchange  information, 
establish  link  */ 

ls=S->L+S->nlink; 

ls->linknr=S->nlink; 

if (commondim<32  &&  ( l<<commondim) <eqmin->nrhs ) 

establishSmallLink (eqmin, eqmax, Is, U+ ( sumnl-commondim) , commondim) 
else 

establishBigLink (eqmin, eqmax, Is , U+ ( sumnl-commondim) , commondim) ; 
S->nlink++; 

} 

for(i=0;  i<sumnl;  ++i) 
free (U [i] ) ; 
free (U) ; 


void  linkSystem ( struct  system  *S) { 

/*  Establishes  connections  for  all  pairs  of  equations  that  can 
exchange  information  */ 
int  i,  j,  eqOni,  eqlni; 
struct  eqSymbol  *eq0,  *eql; 
struct  linkSymbol  *ls; 

S->nlink=0 ; 
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for(i=0;  i<S->neq;  ++i){ 

S->E  [ i ]  . nn=0 ; 

S->E [ i ] . coverID= (u32  ** ) malloc ( ( S->neq-l ) *sizeof (u32  *)); 

S->E [ i ] . link= (struct  links ymbol  ** ) malloc ( ( S->neq-l ) *sizeof (struct 
linkSymbol  *) ) ; 

} 

S->L= (struct  linkSymbol  *) malloc ((( S->neq* ( S->neq- 
1) ) /2) *sizeof (struct  linkSymbol) ) ; 
for(i=0;  i<S->neq-l;  ++i){ 

for(j=i  +  l;  j<S->neq;  ++ j  ) 
tryEstablishLink ( S ,  i ,  j  )  ; 

} 

S->L= (struct  linkSymbol  *) realloc (S->L,  S->nlink*sizeof (struct 
linkSymbol) ) ; 

for(i=0;  i<S->neq;  ++i){ 

S->E [ i] . cover I D= (u32  **)realloc(S->E[i], cover ID, S- 
>E  [i]  . nn*sizeof (u32  *)); 

S->E [i] . link= (struct  linkSymbol  **) realloc (S->E [i] . link, S- 
>E  [i]  .nn*sizeof (struct  linkSymbol  *)); 

} 


void  deLinkSystem ( struct  system  *S) { 

/*  Removes  all  links  in  S  and  frees  the  memory  */ 

int  i; 

struct  eqSymbol  *eq; 

struct  linkSymbol  *li; 

for(eq=S->E;  eq<S->E+S->neq;  ++eq) { 
if (eq->nn>0 ) { 

for(i=0;  i<eq->nn;  ++i) 
free (eq->coverID [ i ] )  ; 
free (eq->coverID) ; 
free (eq->link) ; 
eq->nn=0 ; 

} 

} 

for(li=S->L;  li<S->L+S->nlink;  ++li){ 
free ( li->nCover [ 0 ]  )  ; 
free ( li->nCover [ 1 ] ) ; 
free ( li->cellExists )  ; 

} 

free ( S->L) ; 

S->nlink=0 ; 


void  makeLinkTab ( struct  system  *S,  int  **LT) { 
int  i; 


for(i=0;  i<S->neq;  ++i) 

LT[i]=(int  *)calloc (S->neq, sizeof (int) ) ; 
for(i=0;  i<S->nlink;  ++i){ 

LT [ S->L [ i ] . nlist [ 0 ] - (S->E) ] [S->L[i] . nlist [ 1 ] - ( S- 
>E) ] = (int) ceil (log2 (S->L [i] .ncells) ) ; 
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LT[S->L[i] . nlist [ 1 ] - ( S->E) ] [S->L[i] . nlist [ 0 ] - ( S- 
>E) ]  =  (int) ceil (log2 ( S  —  >L [ i ]  .ncells) ) ; 

} 

} 

void  printLinkTab (int  **LT,  int  n) { 
int  i,  j; 


for(i=0;  i<n;  ++i){ 

for ( j=0;  j<n;  ++j ) { 

printf ("%2d", LT  [i]  [  j  ]  )  ; 
if ( ( j+1) %5==0) 
printf (" | ") ; 

} 

if ( ( i+1 ) %5==0 ) { 
printf ("\n") ; 
for ( j  =0 ;  j  <n;  ++ j ) 
printf  ; 

for(j=0;  j<n/5;  ++ j ) 
printf ; 

} 

printf ("\n")  ; 

} 

} 


83 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


84 


APPENDIX  B 


(These  codes  are  intended  only  to  be  used  for  experiments  with  the 
MRHS  method,  and  neither  their  author  nor  the  author  of  this  thesis  takes  any 
responsibility  for  their  use). 

A.  GENERATION  OF  AN  EQUATION 

/* 

aes_eqs . c 

D.  Canright  2008  Sep  10  Wed  12:44:56 
Generate  MRHS  Equations  for 
Small-scale  Variants  of  the  AES  algorithm 
Also  does  the  encryption! 

Notes:  always  uses  *  form:  last  round  no  MixCols 

always  keysize  =  block  size 
optional  command  line  arguments: 

variant  (string)  =  "nrce"  to  specify  small-scale  variant  of  AES: 


n 

(hex) 

is 

# 

rounds 

(  1  - 

A;  default=A=10  ) 

r 

(int) 

is 

# 

rows  (  1,  2, 

4;  default=4  ) 

c 

(int) 

is 

# 

cols  ( 

CM 

\ — 1 

4;  default=4  ) 

e 

(int) 

is 

# 

bits  in  word 

(  2,  4,  8;  default 

defaults  are  "A448"  for  standard  AES  =  SR* (10;  4;  4; 

8) 

input  (hex)  =  plaintext  block  (default  is  zero  block) 
output  (hex)  =  key  block  (default  is  zero  block) 
outfile  (string)  =  filename  of  output  file  (default  is  stdout) 
while  all  the  above  are  optional,  you  must  have  one  to  have  the  next... 

save  all  the  X  state  data  (output  of  S-box  after  ShiftRows)  and  K  key 
data, 

print  it  out  after  the  equations. 

encryption  re-organized  the  to  give  the  X  state: 

put  ShiftRows  before  S-box,  as  part  of  previous  round 
do  NOT  do  "in  place";  rather,  put  result  in  new  place. 

(ARS)  (MARS)*(n-l)  (A)  rather  than 

(A)  (SRMA)*(n-l)  (SRA)  [where  R  is  RowShift,  S  is  SubstBytes,  ...] 

Does  actual  KeySchedule  and  Encrypt. 

(Note:  keep  InvMix  in  case  do  non-star  versions) . 

*/ 

#include  <stdio.h> 

#include  <string.h> 

#def ine  MAXROUNDS  10 

#def ine  MAXROWS  4 

#def ine  MAXCOLS  4 

#def ine  MAXBITS  8 

#def ine  MAXBLOCK  MAXROWS *MAXC0LS 

#def ine  MAXKEY  MAXBLOCK 

#def ine  MAX VARS  MAXROUNDS *MAXR0WS* (MAXCOLS+1 ) 
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#define  SR(c,r)  (  (c+r)  &  (nCols-1)  ) 


unsigned  char  RoundKeys [ (MAXROUNDS  +1)  *  MAXBLOCK]; 
unsigned  char  States [ (MAXROUNDS  +2)  *  MAXBLOCK]; 
unsigned  int  *Log; 

unsigned  char  *ALog,  *Sbox,  *Mix,  *InvMix,  fieldmask; 

int  nRounds  =  10,  nRows  =  4,  nCols  =  4,  nBits  =  8,  star  =  1,  field, 

block, 

KeyBits,  KeyCols,  nKeyCols; 
int  nEqs,  nVars,  nKeyVars; 

unsigned  char  PT [MAXBLOCK] ,  CT [MAXBLOCK] ,  Eq [2 ] [MAXVARS ] ,  Data[2]; 

enum  InOut  {  In,  Out  }; 

enum  VarType  {  Key,  X,  State  } ; 

unsigned  int  Log8[256]  =  { 


0x00 , 0x00 , 
xDF, 0x03, 

0x19, 

0x01 , 0x32 , 

0x02 , 

OxlA, 

0xC6, 

0x4B, 

0xC7 , 

OxlB, 

0x68, 

0x33, 

OxEE, 

0x64 , 0x04 , 
xlC, OxCl , 

OxEO , 

OxOE, 0x34 , 

0x8D, 

0x81, 

OxEF, 

0x4C, 

0x71, 

0x08 , 

0xC8 , 

0xF8 , 

0x69, 

0x7D, 0xC2 , 
x09, 0x78, 

OxlD, 

0xB5, 0xF9, 

0xB9, 

0x27, 

0x6A, 

0x4D, 

0xE4 , 

0xA6, 

0x72, 

0x9A, 

0xC9, 

0x65 , 0x2F, 
xDA, 0x8E, 

0x8A, 

0x05 , 0x2 1 , 

OxOF, 

OxEl , 

0x24, 

0x12, 

OxFO , 

0x82 , 

0x45, 

0x35 , 

0x93, 

0x96, 0x8F, 
x83, 0x38, 

OxDB, 

OxBD, 0x36, 

OxDO , 

OxCE, 

0x94, 

0x13, 

0x5C, 

0xD2 , 

OxFl , 

0x40, 

0x46, 

0x66, OxDD, 
x91, 0x10, 

OxFD, 

0x30, OxBF, 

0x06, 

0x8B, 

0x62 , 

0xB3 , 

0x25, 

0xE2 , 

0x98, 

0x22 , 

0x88, 

0x7E, 0x6E, 
x3D, OxBA, 

0x48, 

0xC3, 0xA3, 

0xB6, 

OxlE, 

0x42, 

0x3A, 

0x6B, 

0x2  8 , 

0x54, 

OxFA, 

0x85 , 

0x2B, 0x7  9, 
xA7 , 0x57 , 

OxOA, 

0x15 , 0x9B, 

0x9F, 

0x5E, 

OxCA, 

0x4E, 

0xD4 , 

OxAC, 

0xE5 , 

0xF3 , 

0x73, 

OxAF, 0x58, 
xAD, 0xE8 , 

0xA8 , 

0x50 , 0xF4 , 

OxEA, 

0xD6, 

0x74, 

0x4F, 

OxAE, 

0xE9, 

0xD5 , 

0xE7 , 

0xE6, 

0x2C, 0xD7 , 
x51, OxAO, 

0x75, 

0x7A, OxEB, 

0x16, 

OxOB, 

0xF5 , 

0x59, 

OxCB, 

0x5F, 

OxBO , 

0x9C, 

0xA9, 

0x7F, OxOC, 
x7B, 0xB7 , 

0xF6, 

0x6F, 0x17 , 

0xC4 , 

0x49, 

OxEC, 

0xD8 , 

0x43, 

OxlF, 

0x2D, 

0xA4 , 

0x76, 

OxCC, OxBB, 
x29, 0x9D, 

0x3E, 

0x5A, OxFB, 

0x60, 

OxBl , 

0x86, 

0x3B, 

0x52 , 

OxAl , 

0x6C, 

OxAA, 

0x55 , 

0x97 , 0xB2 , 
x5B, OxDl , 

0x87, 

0x90 , 0x61 , 

OxBE, 

OxDC, 

OxFC, 

OxBC, 

0x95, 

OxCF, 

OxCD, 

0x37, 

0x3F, 

0x53, 0x39, 
xD3 , OxAB, 

0x84, 

0x3C, 0x41, 

0xA2 , 

0x6D, 

0x47, 

0x14, 

0x2A, 

0x9E, 

0x5D, 

0x56, 

0xF2 , 

0x44, 0x11, 
xE3 , 0xA5 , 

0x92, 

0xD9, 0x23, 

0x20, 

0x2E, 

0x89, 

0xB4 , 

0x7C, 

0xB8 , 

0x26, 

0x77, 

0x99, 

0x67 , 0x4A, 
x7  0 , 0x07 , 
}; 

OxED, 

OxDE, 0xC5 , 

0x31, 

OxFE, 

0x18, 

OxOD, 

0x63, 

0x8C, 

0x80, 

OxCO , 

0xF7 , 

unsigned  char  ALog8 [256] 

=  ( 

0x01, 0x03, 
xl3, 0x35, 

0x05 , 

OxOF, Oxll , 

0x33, 

0x55 , 

OxFF, 

OxlA, 

0x2E, 

0x72, 

0x96, 

OxAl, 

0xF8 , 

0x5F, OxEl , 

0x38 , 

0x4  8 , 0xD8 , 

0x73, 

0x95, 

0xA4 , 

0xF7 , 

0x02 , 

0x06, 

OxOA, 

OxlE, 

0x22 , 

x66, OxAA, 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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0xE5 , 0x34 , 0x5C, 
xE6, 0x31, 

0xE4 , 

0x37, 

0x59, OxEB, 

0x26, 

0x6A, 

OxBE, 0xD9 , 

0x70, 

0x90 

0x53 , 0xF5 , 0x04  , 
xB2 , OxCD, 

OxOC, 

0x14, 

0x3C, 0x44 , 

OxCC, 

0x4F, 

OxDl , 0x68 , 

0xB8 , 

0xD3 

0x4C, 0xD4, 0x67, 
x78, 0x88, 

0xA9, 

OxEO , 

0x3B, 0x4D, 

0xD7 , 

0x62 , 

0xA6, OxFl , 

0x08 , 

0x18 

0x83, 0x9E, OxB9, 
x76, 0x9A, 

OxDO , 

0x6B, 

OxBD, OxDC, 

0x7F, 

0x81, 

0x98 , 0xB3 , 

OxCE, 

0x49 

0xB5, 0xC4, 0x57, 
x61 , 0xA3 , 

0xF9, 

0x10, 

0x30 , 0x50 , 

OxFO, 

OxOB, 

OxlD, 0x27 , 

0x69, 

OxBB 

OxFE, 0x19, 0x2B, 
x60 , OxAO , 

0x7D, 

0x87, 

0x92 , OxAD, 

OxEC, 

0x2F, 

0x7 1 , 0x93 , 

OxAE, 

OxE9 

OxFB, 0x16, 0x3A, 
x3F, 0x4 1 , 

0x4E, 

0xD2 , 

0x6D, 0xB7 , 

0xC2 , 

0x5D, 

0xE7 , 0x32 , 

0x56, 

OxFA 

0xC3, 0x5E, 0xE2, 
xDA, 0x7  5 , 

0x3D, 

0x47, 

0xC9, 0x40, 

OxCO, 

0x5B, 

OxED, 0x2C, 

0x74, 

0x9C 

0x9F, OxBA, 0xD5 , 
x89, 0x80, 

0x64, 

OxAC, 

OxEF, 0x2A, 

0x7E, 

0x82 , 

0x9D, OxBC, 

OxDF, 

0x7A 

0x9B, 0xB6, OxCl , 
xC5 , 0x54 , 

0x58 , 

0xE8 , 

0x23, 0x65, 

OxAF, 

OxEA, 

0x25, 0x6F, 

OxBl , 

0xC8 

OxFC, OxlF, 0x21, 
x46, OxCA, 

0x63, 

0xA5 , 

0xF4 , 0x07 , 

0x09, 

OxlB, 

0x2D, 0x77 , 

0x99, 

OxBO 

0x45, OxCF, 0x4A, 
xF3 , OxOE, 

OxDE, 

0x79, 

0x8B, 0x86, 

0x91, 

OxA8 , 

0xE3, 0x3E, 

0x42, 

0xC6 

0x12, 0x36, 0x5A, 
xOD, 0x17 , 

OxEE, 

0x29, 

0x7B, Ox8D, 

0x8C, 

Ox8F, 

0x8A, 0x85, 

0x94, 

0xA7 

0x39, 0x4B, OxDD, 
xF6, 0x01, 

0x7C, 

0x84, 

0x97 , 0xA2 , 

OxFD, 

OxlC, 

0x24 , 0x6C, 

0xB4 , 

0xC7 

}  ; 

unsigned  int  Log4[16]  =  { 

0,  0,  1,  4,  2,  8,  5,  10,  3,  14,  9,  7,  6,  13,  11,  12, 

}  ; 

unsigned  char  ALog4[16]  =  { 

1,  2,  4,  8,  3,  6,  12,  11,  5,  10,  7,  14,  15,  13,  9,  1, 

}  ; 

unsigned  int  Log2[4]  =  { 

0,  0,  1,  2, 

}  ; 

unsigned  char  ALog2[4]  =  { 

1,  2,  3,  1, 

}  ; 

unsigned  char  Sbox8 [256]  =  { 

0x63, 0x7C,  0x77, 0x7B, 0xF2, 0x6B, 0x6F, 0xC5, 0x30, 0x01, 0x67, 0x2B, OxFE, 
xAB, 0x7  6, 

OxCA, 0x82, 0xC9, 0x7D, OxFA, 0x59, 0x47, OxFO, OxAD, 0xD4, 0xA2, OxAF, 0x9C, 
x72, OxCO, 

0xB7,  OxFD,  0x93, 0x26,  0x36,  0x3F, 0xF7, OxCC, 0x34, 0xA5, 0xE5, OxFl, 0x71, 
x31, 0x15, 

0x04, 0xC7, 0x23, 0xC3, 0x18, 0x96,  0x05, 0x9A,  0x07, 0x12, 0x80, 0xE2, OxEB, 
xB2 , 0x7  5 , 


OxAB, 0 
0x6E, 0 
0x2  8 , 0 
OxDB, 0 
0xD6, 0 
0x20, 0 
0x15,0 
OxBF, 0 
0x8E, 0 
0x43,0 
OxCB, 0 
0x51,0 
0xF2 , 0 
0x52 , 0 


0xD7 , 0 
0xA4 , 0 
0xD8 , 0 
0x27,0 
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0x09, 0x83, 0x2C, OxlA, OxlB, 0x6E, 0x5A, OxAO, 0x52, 0x3B, 0xD6,  0xB3,  0x29, 0xE3, 0 
x2F, 0x84 , 

0x53, OxDl, 0x00, OxED, 0x20, OxFC, OxBl, 0x5B, 0x6A, OxCB, OxBE, 0x39, 0x4A, 0x4C, 0 
x58, OxCF, 

OxDO, OxEF, OxAA, OxFB, 0x43, 0x4D, 0x33, 0x85, 0x45, 0xF9, 0x02, 0x7F, 0x50, 0x3C, 0 
x9F, 0xA8 , 

0x51, 0xA3, 0x40, 0x8F, 0x92, 0x9D, 0x38, 0xF5, OxBC, 0xB6, OxDA, 0x21, 0x10, OxFF, 0 
xF3 , 0xD2 , 

OxCD, OxOC, 0x13, OxEC, 0x5F,  0x97, 0x44, 0x17, 0xC4, 0xA7, 0x7E,  0x3D,  0x64, 0x5D,  0 
xl9, 0x73, 

0x60, 0x81, 0x4F, OxDC, 0x22, 0x2 A,  0x90, 0x88, 0x46, OxEE, 0xB8, 0x14, OxDE,  0x5E,  0 
xOB, OxDB, 

OxEO,  0x32, 0x3A, OxOA,  0x49,  0x06,  0x24, 0x5C, 0xC2, 0xD3, OxAC,  0x62, 0x91, 0x95, 0 
xE4 , 0x7  9 , 

0xE7,  0xC8, 0x37, 0x6D,  Ox8D,  0xD5, 0x4E, 0xA9, 0x6C, 0x56, 0xF4, OxEA, 0x65, 0x7A, 0 
xAE, 0x08, 

OxBA,  0x78, 0x25, 0x2E,  OxlC,  0xA6,  0xB4, 0xC6,  OxE8,  OxDD,  0x74,  OxlF,  0x4B,  OxBD,  0 
x8B, Ox8A, 

0x7  0, 0x3E, 0xB5, 0x66,  0x48, 0x03, 0xF6, OxOE,  0x61, 0x35, 0x57, 0xB9,  0x86,  OxCl,  0 
xlD, 0x9E, 

OxEl, 0xF8, 0x98, Oxll, 0x69, 0xD9, Ox8E, 0x94, 0x9B, OxlE, 0x87, 0xE9, OxCE, 0x55, 0 
x2  8 , OxDF, 

0x8C, OxAl, 0x89, OxOD, OxBF, 0xE6, 0x42, 0x68, 0x41, 0x99, 0x2D, OxOF, OxBO, 0x54, 0 
xBB, 0x16, 

}  ; 

unsigned  char  Sbox4[16]  =  { 

0x6,  OxB,  0x5, 0x4, 0x2,  OxE, 0x7, OxA, 0x9, OxD, OxF, OxC, 0x3, Oxl, 0x0, 0x8, 

}  ; 


unsigned  char  Sbox2[4]  =  { 

2,  3,  1,  0, 

}  ; 


unsigned  char  Mix4[4]  =  { 
0x2 , 0x3, Oxl, Oxl, 

}  ; 


unsigned  char  InvMix4[4]  =  { 
OxE, OxB, OxD, 0x9, 

}  ; 

unsigned  char  InvMix42[4]  =  { 
0x0, 0x2, 0x3, 0x0, 

}  ; 

unsigned  char  Mix2[2]  =  { 

0x3 , 0x2 , 

}  ; 

unsigned  char  Mixl[l]  =  { 

Oxl, 

}  ; 


//  multiply  by  "2"  in  field 
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#define  P0LY8  OxlB 
#define  POLY4  0x13 
#define  P0LY2  0x07 
unsigned  char  HIBIT,  POLY; 
unsigned  char  mul2  (unsigned  char  x)  { 
unsigned  char  y; 
y  =  x  <<  1 ; 

if  (  x  &  HIBIT  )  y  A=  POLY; 
return (  y  ) ; 

} 

//  multiply  two  bytes  in  field 

unsigned  char  mul (unsigned  char  x,  unsigned  char  y) 

{ 

if  (x  &&  y) 

return  (ALog [ (Log [x]  +  Log[y])  %  (field  -  1)]); 

else 

return  (  0 ) ; 

} 

#include  "eqs__io.h"  //  include  I/O  package 
//  set  up  specific  small-scale  variant  of  AES 

//  assumes  main ()  already  set:  nRounds,  nRows,  nCols,  nBits,  star 
int  setup (void) 

{ 

int  returnval  =  0; 

//  check  parameters  for  validity 

if  (nRounds  <  1  | |  nRounds  >  MAXROUNDS)  { 
nRounds  =  10; 
returnval  =  1; 

} 

if  (  !  (nCols  ==  1  |  |  nCols  ==  2  |  |  nCols  ==  4 ) )  { 

nCols  =  4; 
returnval  =  1; 

} 

switch  (nBits)  { 
case  2 : 

Log  =  Log2; 

ALog  =  ALog2; 

Sbox  =  Sbox2; 

POLY  =  POLY2 ; 
field  =  4; 
break; 
case  4  : 

Log  =  Log4; 

ALog  =  ALog4; 

Sbox  =  Sbox4; 

POLY  =  POLY4 ; 
field  =  16; 
break; 
default : 
nBits  =  8; 
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returnval  =1;  //  if  bad  value,  use  default,  fall  thru 

case  8  : 

Log  =  Log8; 

ALog  =  ALog8; 

Sbox  =  Sbox8; 

POLY  =  POLY 8 ; 
field  =  256; 
break; 

} 

switch  (nRows)  { 
case  1 : 

Mix  =  InvMix  =  Mixl; 
break; 
case  2  : 

Mix  =  InvMix  =  Mix2; 
break; 
default : 
nRows  =  4 ; 

returnval  =1;  //if  bad  value,  use  default,  fall  thru 

case  4  : 

Mix  =  Mix4 ; 

InvMix  =  (nBits  ==  2)  ?  InvMix42  :  InvMix4; 

break; 

} 

fieldmask  =  field  -  1; 

HIBIT  =  1  «  (nBits  -  1); 

setScale();  //  set  up  bit  matrices  for  scalars 
block  =  nRows  *  nCols; 

KeyBits  =  block  *  nBits; 
nKeyCols  =  (nRounds  +1)  *  nCols; 
nKeyVars  =  block  +  nRounds  *  nRows; 
nVars  =  nKeyVars  +  block  *  (nRounds  -  1); 
nEqs  =  nVars; 

return  returnval; 

} 

int  KeySchedule (unsigned  char  Key[]) 

{ 

int  colbits,  returnval  =  0; 
int  r,  c; 

unsigned  char  col [MAXROWS ] ,  t,  rcon; 

colbits  =  nRows  *  nBits; 

KeyCols  =  KeyBits  /  colbits; 

#def ine  NOISY  0 
#if  NOISY 

fprintf (stderr, "-KeySched:  colbits=%d,  KeyCols=%d,  nKeyCols=%d, 
Key=%p\n" , 

colbits,  KeyCols,  nKeyCols,  Key) ; 
fprintf (stderr, "-Key :  "); 

for  (r  =0;  r  <  block;  r++)  fprintf ( stderr ,  ((nBits>4)  ?  "%02X"  : 
"%01X"),  Key [r] ) ; 

fprintf (stderr, "\n" ) ;  fflush (stderr) ; 
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#endif 

/*  Copy  key  */ 

for  (c  =0;  c  <  KeyCols;  C++) 
for  (r  =  0;  r  <  nRows;  r++) 

RoundKeys[r  +  nRows  *  c]  =  Key[r  +  nRows  *  c] ; 
for  (r  =  0;  r  <  nRows;  r++) 

col[r]  =  Key[r  +  nRows  *  ( c  —  1 ) ] ; 

#if  NOISY 

fprintf (stderr, "-KeySched:  Key  copied;  c=%d,  col=  ",  c)  ; 

for  (r  =0;  r  <  nRows;  r++)  fprintf ( stderr ,  (nBits>4 ) ?"%02X" : "%01X" , 
col [r] ) ; 

fprintf (stderr, "\n" ) ;  fflush (stderr) ; 

#endif 

for  (rcon  =1;  c  <  nKeyCols;  C++)  { 

/*  calculate  new  columns  until  enough  */ 
if  (c  %  KeyCols  ==  0)  { 

t  =  col [  0  ]  ; 

for  (r  =0;  r  <  (nRows  -  1) ;  r++) 
col [r]  =  Sbox[col[r  +  1 ] ] ; 
col [nRows  -  1]  =  Sbox[t]; 
col [ 0 ]  A=  rcon; 

#if  NOISY 

fprintf (stderr, "-KeySched:  apply  F;  t=%X,  rcon=%X,  col=  ",  t,  rcon); 

for  (r  =0;  r  <  nRows;  r++)  fprintf ( stderr ,  (nBits>4 ) ?"%02X" : "%01X" , 
col [r] ) ; 

fprintf (stderr, "\n" ) ;  fflush (stderr) ; 

#endif 

rcon  =  mul(2,  rcon); 

} 

//  need  to  handle  KeyCols  =  1  differently 
for  (r  =  0;  r  <  nRows;  r++) 

RoundKeys[r  +  nRows  *  c]  =  (KeyCols  ==  1)  ?  col [r]  : 

(col [r]  A=  RoundKeys[r  +  nRows  *  (c  -  KeyCols)]); 

#if  0 

fprintf (stderr, "-KeySched:  RoundKeys  col[%d]=  ",  c)  ; 
for  (r  =0;  r  <  nRows;  r++)  fprintf ( stderr ,  (nBits>4 ) ?"%02X" : "%01X" , 
RoundKeys [r  +  nRows  *  c] )  ; 
fprintf (stderr, " \n" ) ;  fflush (stderr) ; 

#endif 

} 

return  returnval; 

} 

//  do  one  round  on  block:  (ARS)  for  #0  or  else  (MARS) 
void  doround (unsigned  char  State [],  unsigned  char  roundKey[], 
int  round) 

{ 

unsigned  char  t[MAXROWS]; 
int  i,  r,  c,  offset=0; 

if  (round)  //  if  normal  round 

for  (c  =0;  c  <  nCols;  C++)  { 

for  (r  =  0;  r  <  nRows;  r++) 

for  (t[r]  =  i  =  0;  i  <  nRows;  i++) 
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t[r]  A  =  //  MixColumns 

mul(Mix[i],  State [ ( (r  +  i)  %  nRows)  +  nRows  * 

c]  )  ; 

for  (r  =  0;  r  <  nRows;  r++) 

State [block  +  r  +  nRows  *  c]  =  t[r]; 

} 

else  offset  =  -block; 

State  +=  block; 

for  (i  =0;  i  <  block;  i++) 

State [i]  =  State [offset+i]  A  roundKey[i];  //  AddRoundKey 

for  (r  =  1;  r  <  nRows;  r++)  { 

for  (c  =0;  c  <  nCols;  C++)  //  ShiftRows 

t[c]  =  State [r  +  nRows  *  ( (c  +  r)  %  nCols)]; 

for  (c  =0;  c  <  nCols;  C++) 

State [r  +  nRows  *  c]  =  t[c]; 

} 

for  (i  =0;  i  <  block;  i++) 

State [i]  =  Sbox [ State  [ i ]] ;  //  SubBytes 


} 

//  do  round  #n  on  block:  (A) 

void  doroundn (unsigned  char  State [],  unsigned  char  roundKey[]) 

{ 

int  i; 

for  (i  =0;  i  <  block;  i++) 

State [block  +  i]  =  State [i]  A  roundKey[i];  // 

AddRoundKey 

} 

//  encrypt  block  (NOT  in  place  -  keep  output  of  each  S-box) 
void  encrypt (  void  ) 

{ 

int  i,  round; 

for  (i  =0;  i  <  block;  i++) 

States [i]  =  PT[i];  //  copy  PT  in 

for  (round  =  0;  round  <  nRounds;  round++)  { 

doround (  States  +  round*block,  RoundKeys  +  round*block,  round); 

} 

doroundn ( States  +  round*block,  RoundKeys  +  round*block) ; 
for  (i  =  0;  i  <  block;  i++) 

CT[i]  =  States [ (nRounds+1 ) *block  +  i];  //  copy  CT  out 

} 

void  NewEq (  void  ) 

{ 

int  i,  r; 

for  (r  =  In;  r  <=  Out;  r++)  { 

Data[r]  =  0; 

for  (i  =  0;  i  <  nVars;  i++) 

Eq [ r ] [i]  =  0; 

} 
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} 


int  VarNum (  enum  VarType  var, 

int  round,  int  col,  int  row  ) 

{ 

switch  (var)  { 
case  Key: 

if  (round  ==  0) 

return  (  col  *  nRows  +  row  ) ; 
else  / /  then  col  ==  0 

return  (  block  +  (round-1)  *  nRows  +  row  ); 
case  X: 

return  (  nKeyVars  +  (round-1)  *  block  +  col  *  nRows  +  row  ) ; 

} 

return  (  0  ) ;  //  dummy 

} 

void  AddVar (  enum  InOut  line,  enum  VarType  var, 
int  round,  int  col,  int  row,  int  scale  ) 

{ 

int  r; 

switch  (var)  { 
case  Key: 

if  (round  ==  0  |  col  ==  0) 

Eq[line] [  VarNum (  Key,  round,  col,  row  )  ]  A=  scale; 

else  { 

AddVar (  line.  Key,  round,  col-1,  row,  scale  ); 

AddVar (  line.  Key,  round-1,  col,  row,  scale  ); 

} 

break; 
case  X: 

Eq[line] [  VarNum (  X,  round,  col,  row  )  ]  A=  scale; 

break; 

case  State:  //  scale  must  be  1 
for  (r  =  0;  r  <  nRows;  r++)  { 

AddVar (  line,  X,  round,  col,  (row+r) & (nRows-1) ,  Mix[r]  ) 

} 

AddVar (  line.  Key,  round,  col,  row,  1  ); 
break; 

} 

} 

void  KeyScheduleEqs (void) 

{ 

int  i,  r; 

unsigned  char  rcon; 

for  (i  =  1,  rcon  =1;  i  <=  nRounds;  i++)  { 

for  (r  =  0;  r  <  nRows;  r++)  { 

NewEq ( ) ; 

AddVar (  Out,  Key,  i,  0,  r,  1); 

if  (  nCols  >  1  )  AddVar (  Out,  Key,  i-1,  0,  r,  1); 

AddVar (  In,  Key,  i-1,  nCols-1,  (r+1) & (nRows-1) ,  1); 
if  (  r  ==  0  )  Data [  Out  ]  =  rcon; 
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WriteEq ( ) ; 


} 

rcon  =  mul2 (rcon) ; 

} 

} 

//  do  only  1  round  on  block 
void  doonlyroundEqs ( int  round) 

{ 

int  r,  c; 

for  (c  =0;  c  <  nCols;  C++)  { 

for  (r  =  0;  r  <  nRows;  r++)  { 

NewEq ( ) ; 

AddVar (  Out,  Key,  round,  c,  r,  1); 

Data [  Out  ]  =  CT [  r  +  nRows  *  c  ] ; 
AddVar (  In,  Key,  round-1,  SR(c,r),  r,  1) 
Data [  In  ]  =  PT [  r  +  nRows  *  SR(c,r)  ]; 
WriteEq ( ) ; 

} 

} 

} 

//  do  round  #1  on  block 
void  doroundlEqs ( int  round) 

{ 

int  r,  c; 

for  (c  =0;  c  <  nCols;  C++)  { 

for  (r  =  0;  r  <  nRows;  r++)  { 

NewEq ( ) ; 

AddVar (  Out,  X,  round,  c,  r,  1); 

AddVar (  In,  Key,  round-1,  SR(c,r),  r,  1) 
Data [  In  ]  =  PT [  r  +  nRows  *  SR(c,r)  ]; 
WriteEq ( ) ; 

} 

} 

} 

//  do  one  round  on  block 
void  doroundEqs ( int  round) 

{ 

int  r,  c; 

for  (c  =0;  c  <  nCols;  C++)  { 

for  (r  =  0;  r  <  nRows;  r++)  { 

NewEq ( ) ; 

AddVar (  Out,  X,  round,  c,  r,  1); 

AddVar (  In,  State,  round-1,  SR(c,r),  r, 
WriteEq ( ) ; 

} 

} 

} 

//  do  round  #n  on  block 
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void  doroundnEqs ( int  round) 

{ 

int  r,  c; 

for  (c  =0;  c  <  nCols;  C++)  { 

for  (r  =  0;  r  <  nRows;  r++)  { 

NewEq ( ) ; 

AddVar (  Out,  Key,  round,  c,  r,  1); 

AddVar (  In,  State,  round-1,  SR(c,r),  r,  1); 
Data [  Out  ]  =  CT [  r  +  nRows  *  c  ]  ; 

WriteEq ( ) ; 

} 

} 

} 

void  EncryptEqs (void) 

{ 

int  round; 


if  (  nRounds  ==  1  )  { 

doonlyroundEqs (1) ; 
return; 

} 

doroundlEqs (1) ; 

for  (round  =2;  round  <  nRounds;  round++)  { 
doroundEqs (round) ; 

} 

doroundnEqs (round) ; 


int  main (int  argc,  char  *argv[]) 

{ 

unsigned  char  Key [MAXBLOCK] ; 


if  (argc  >  1)  { 

sscanf (argv[l] ,  "%lx%ld%ld%ld" , 

SnRounds,  SnRows,  SnCols,  SnBits); 

} 

fprintf ( stderr ,  "  nRounds=%d,  nRows=%d,  nCols=%d,  nBits 
star=%d\n" , 

nRounds,  nRows,  nCols,  nBits,  star); 
if  ( setup ( ) ) 

fprintf (stderr, 

"Bad  parameter (s) ;  now:\n  nRounds=%d,  nRows=%d, 
nBits=%d,  star=%d\n", 

nRounds,  nRows,  nCols,  nBits,  star) ; 

//  by  default  KeyBits  =  bits  in  block 

ReadBlock(  (argc  >  2)  ?  argv[2]  :  PT) ; 

ReadBlock(  (argc  >  3)  ?  argv[3]  :  Key)  ; 

if  (argc  >  4) 

if  ( f reopen (argv [ 4 ] ,  "w",  stdout)  !=  stdout)  { 

fprintf ( stderr ,  "Could  not  open  output  file  %s\n" 
return  1; 


} 


%d. 


nCols=%d, 


argv [4] ) 
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KeySchedule (Key) ; 
encrypt ( ) ; 

WriteSystemHeader ( ) ; 

KeyScheduleEqs () ; 

//  do  only  single  block 
EncryptEqs ( ) ; 

WriteVars ( ) ; 

printf("  nRounds=%d,  nRows=%d,  nCols=%d,  nBits=%d,  star=%d\n", 
nRounds,  nRows,  nCols,  nBits,  star); 

WriteKeys ( )  ; 

WriteStates  ( )  ; 

return  (0) ; 

} 

B.  MRHS  algorithm 

/* 

mrhs . c 

D.  Canright  2008  Aug  22  Fri  13:58:27 
Interactive  tool  to  solve  MRHS  equations 

input  file  format: 
header : 

NVAR  #  variables 
neq  #  equations  (symbols) 
for  each  equation  symbol: 
nlin  #  rows 
nrhs  #  RHS 

A  (by  rows,  binary) 

b  (by  cols,  binary) 

*/ 

#include  "mrhs.h"  //  includes  all  others 
#include  <ctype.h> 

int  main(int  argc,  char  *argv[] ) 

{ 

int  i,  n; 

char  input [101],  *filename; 
char  menu [ ]  =  { 

"Enter  commands  from  the  menu  below. \n" 

"  #  means  a  number  is  required:  use  nondigit  for  default  value  (or  for 

'all' ) \n" 

"1  linkSystem ( S ) ; \n" 

"a  agreeSystem ( S ) ; \n" 

"d  deLinkSystem ( S ) ; \n" 

"g  #  packSystem ( S ,  #) ;  (glue) \n" 

"x  extractLinearlnf o ( S ) ; \n" 

"w  writeSystem ( S ) ; \n" 

"i  #  Equationlnf o ( &S->E [ # ] ) ; \n" 
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"p  #  printEquation ( &S->E [ # ] ) ; \n" 

"e  printLinEquation ( S->linbank,  0,  S->linbank->nlin) ; \n" 

"t  writeLinkTab ( S ,  36);\n" 

"q  QUIT  (exit) \n" 

}  ; 

//  struct  eqSymbol  *tmpeq  =  (struct  eqSymbol  *)  malloc (sizeof (struct 
eqSymbol) ) ; 

struct  system  *S  =  (struct  system  *)  malloc (sizeof (struct  system 

) )  ; 

CHECKPTR ( S ) ;  //  test  of  macro 

/* 

read  in  header,  set  up  sys,  setup  &  read  eqns 
link  system 
agree 
glue 
*/ 

#if  0 

i  =  0; 

CHECKPTR (i);  //  test  of  macro 

return ( 0 )  ; 

#endif 

if  (  argc  >  1  )  filename  =  argv[l]; 
else  { 

fprintf (  stderr,  "Enter  filename  for  input  system:  "); 
scanf ( "%100s" ,  input); 
filename  =  input; 

} 

if  (  readSystem ( S ,  filename)  )  { 

fprintf (  stderr,  "Error  in  main:  did  not  get  input  system;  abort !\n" 

)  ; 

return  (  1  )  ; 

} 

printf ( "got :  NVAR  =  %d  (NWORDS  =  %d) ;  neq  =  %d\n",  NVAR,  NWORDS,  S->neq 

)  ; 

printf (menu) ; 

while  (  scanf ( "%100s" , input)  ==  1  )  { 

switch  (input  [0])  { 

case  ' a ' : 
case  ' A ' : 

printf  ("**  agreeSystem\n" ) ; 

agreeSystem ( S ) ; 

break; 

case  ' e ' : 

case  ' E ' : 

printLinEquation ( S->linbank,  0,  S->linbank->nlin) ; 

break; 

case  ' d ' : 

case  ' D ' : 

printf ("**  deLinkSystem\n"); 
deLinkSystem ( S ) ; 
break; 
case  ' p ' : 
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case  ' P ' : 

if  (  input  [1]  )  i=l; 

else  {  scanf ( "%100s" , input) ;  i  =  0;  } 

if  (  isdigit ( input  [ i ] )  )  { 

sscanf (  input+i,  "%d",  &i  ); 

if  (  i  <  0  ||  i  >=  S->neq  )  { fprintf ( stderr , "  bad  eq  #:  %d\n",i); 

break; } 

printEquation ( &S->E [ i ]  )  ; 

} 

else 

for(i=0;  i<S->neq;  ++i)  printEquation ( &S->E [ i ]) ; 
break; 
case  ' i  '  : 
case  '  I '  : 

if  (  input  [1]  )  i=l; 

else  {  scanf ( "%100s ",  input ) ;  i  =  0;  } 

if  (  isdigit ( input  [ i ] )  )  { 

sscanf (  input+i,  "%d",  &i  ); 

if  (  i  <  0  ||  i  >=  S->neq  )  { fprintf ( stderr , "  bad  eq  #:  %d\n",i); 

break; } 

Equation Info ( &S->E [ i ]  )  ; 

} 

else 

for(i=0;  i<S->neq;  ++i)  Equationlnf o ( &S->E [ i ] ) ; 
break; 
case  ' g ' : 
case  ' G ' : 

if  (  input  [1]  )  i=l; 

else  {  scanf ( "%100s" , input) ;  i  =  0;  } 

if  (  isdigit ( input  [ i ] )  ) 

sscanf (  input+i,  "%d",  &i  ); 
else 

i  =  1024;  //  arbitrarily  picked  a  (small)  max  size  of  glued  eqns 

printf("**  packSystem  (glue);  max  =  %d\n",  i) ; 
packSystem ( S ,  i) ; 
break; 
case  ' 1 '  : 
case  ' L ' : 

printf("**  linkSystem\n" ) ; 

linkSystem ( S ) ; 

break; 

case  ' q ' : 

case  ' Q ' : 

return  0; 

case  ' t ' ; 

case  ' T ' : 

writeLinkTab  ( S ,  36); 

break; 

case  ' w ' : 

case  ' W ' : 

writeSystem ( S ) ; 

break; 

case  ' x ' : 

case  ' X ' : 

printf ("**  extractLinearInfo\n" ) ; 
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extractLinearlnf o (S) ; 

break; 

default : 

printf (menu) ; 

break; 

} 

printf ("now:  neq  =  %d  ;  nlink  =  %d  ;  linbank  =  %d\n", 
S->neq,  S->nlink,  S->linbank->nlin  ) ; 

} 

return  (0) ; 

} 


99 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


100 


LIST  OF  REFERENCES 


[1]  W.  Trappe  and  L.C.  Washington,  Introduction  to  Cryptography  with  Coding 
Theory,  2d  ed.,  Upper  Saddle  River:  Pearson  Prentice  Hall,  2006. 

[2]  J.  Daemen  and  V.  Rijmen,  The  Design  of  Rijndael:  AES — The  Advanced 
Encryption  Standard,  New  York:  Springer,  2002. 

[3]  “Announcing  the  Advanced  Encryption  Standard,”  Federal  Information 
Processing  Standards  Publication  197,  November  26,  2001,  [Online], 
Available:  http://csrc.nist.gov/publications/fips/fips1 97/fips-1 97.pdf 
(accessed  June  2008). 

[4]  C.  Cid,  S.  Murphy  and  M.  Robshaw,  Algebraic  Aspects  of  the  Advanced 
Encryption  Standard,  New  York:  Springer,  2006. 

[5]  H.  Raddum  and  I.  Semaev,  “Solving  MRHS  linear  equations,”  Cryptology 
ePrint  Archive,  Report  2007/285,  August  28,  2007, 
http://eprint.iacr.org/2007/285  (accessed  September  2008). 

[6]  H.  Raddum  and  I.  Semaev  (private  communication:  provided  codes  in 
Appendix  A). 

[7]  “Advanced  encryption  standard,”  Wikipedia,  September  2008.  [Online], 
Available:  http://en.wikipedia.org/wiki/Advanced  Encryption  Standard 
(accessed  September  2008). 

[8]  E.  Trichina,  L.  Korkishko,  Secure  and  Efficient  AES  Software 
Implementation  for  Smart  Cards,  [Online],  Available: 
http://eprint.iacr.org/2004/149.pdf  (accessed  June  2009). 

[9]  C.  Cid,  S.  Murphy  and  M.J.B.  Robshaw,  “Small-scale  variants  of  the  AES,” 
Fast  Software  Encryption:  12th  International  Workshop  (FSE  2005),  Paris, 
France,  February  21-23,  2005,  Lecture  Notes  in  Computer  Science,  vol. 
3557,  pp.  145-162,  Berlin:  Springer,  2005. 


101 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


102 


INITIAL  DISTRIBUTION  LIST 


1.  Defense  Technical  Information  Center 
Ft.  Belvoir,  Virginia 

2.  Dudley  Knox  Library 
Naval  Postgraduate  School 
Monterey,  California 

3.  Professor  David  Canright 
Department  of  Applied  Mathematics 
Naval  Postgraduate  School 
Monterey,  California 


103 


